diff --git a/CHANGELOG.md b/CHANGELOG.md index 03fd58d..bf3eb30 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,11 +20,40 @@ and this project adheres to [Semantic Versioning][]. `tl.external.scores_to_gmt`/`scores_to_covar`, `tl.external.compute_escore`, `JointNMFWrapper`) - New "DonorData basics" tutorial: building a `DonorData` from your own genotype/expression data (including loading genotypes via `read_sgkit_zarr`/`stream_pgen_to_zarr`), how donor - syncing works, selecting subsets, aggregating, and saving/loading — no analysis, just the + syncing works, selecting subsets, aggregating, and saving/loading. No analysis, just the data structure itself +- `tl.external.build_known_cis_eqtls_from_tensorqtl`: build a known-cis-eQTL annotation + (variant x gene, binary) from a completed TensorQTL nominal cis-scan, for use as a + fine-mapping prior ### Fixed +- `DonorData.copy()` built a genuinely new object but always copied `_G`/`_C` + regardless of whether they were views, unlike its previous behavior; reverted to + only copying `_G`/`_C` when they're actually views (mutating `self` and returning + it, as before) +- `io.read_zarr_dd` eagerly materialized `G`/`C` on load, unlike `io.read_pgen_zarr` + (which already kept a dense `X` Dask-backed); it now uses the same lazy read path, + so selecting a subset (e.g. one cell type of `C`) never touches the rest of the data +- `DonorData.write_zarr_dd` picked its own chunk shape for a dense `X` without regard + to any chunking the input already had (confirmed: it disregarded a Dask array's own + `.chunks` entirely), which for a genome-scale `X` could pick a shape badly misaligned + with the data's on-disk layout and make the write far slower than necessary; `X` is + now written separately with sane chunking (a Dask array's own chunks, or capped at + 4096 per axis for a plain array), overridable via `write_zarr_dd`'s new `x_chunks` + argument; the same fix now also applies to any dense `layers` entry (previously only + `X` was chunked sanely, so a dense layer the same shape as `X` (e.g. raw counts kept + alongside a normalized `X`) still got the pathological chunking) +- `io.read_zarr_dd`'s lazy read only kept `X` Dask-backed, not a dense `layers` entry, + so a dense layer was still eagerly materialized in full on load (and could exhaust + memory for a genome-scale one); `layers` is now read lazily the same way `X` is +- `io._pgen._read_pvar` only mapped a `.pvar`'s columns to cellink's canonical + `chrom`/`pos`/`snp_id`/`a0`/`a1` names when a `#CHROM` header line was present; + a headerless `.pvar` now gets the same canonical mapping, assuming the standard + PLINK2 column order, instead of being left with unnamed columns +- `io.from_sgkit_dataset` didn't cast `chrom`/`a0`/`a1` to `str`, unlike the pgen + reader, which could leave a numeric `chrom` dtype that breaks a downstream string + membership check (e.g. in tensorQTL) - `pp.log_transform` and `pp.normalize` no longer mutate the input `AnnData` when called with `inplace=False` - `Skat` and `StructLMM` were unusable (missing `limix-core` imports); both now work, @@ -41,6 +70,24 @@ and this project adheres to [Semantic Versioning][]. contributing guide - A stray unanchored `data` entry in `.gitignore` silently hid any new file added under `tests/data/` or `docs/tutorials/data/`; scoped it to the latter only +- `io.to_plink` crashed (`TypeError: NDFrame.to_csv() got an unexpected keyword + argument 'line_terminator'`) with `pandas-plink` 2.2.9 on a modern pandas, which + renamed/removed that `to_csv` kwarg; bumped the minimum `pandas-plink` version to + 2.3.0, which uses the current `lineterminator` kwarg internally + +### Known limitations + +- Reading a sparse (`csr`/`csc`) `X` or layer back from a Zarr v3 store is slow: + confirmed via two independent hand-built-matrix tests that plain + `anndata.io.read_elem` alone (no cellink code involved) takes ~19s for a + 5M-nnz array that writes in ~0.2s, and `DonorData.write_zarr_dd`/`read_zarr_dd` + add only ~9% on top of that baseline, so this is inherent to AnnData's + sparse Zarr v3 read path in this environment, not a `DonorData`-specific + issue. Not yet root-caused further, and not currently fixed, since no real + data in this codebase's own pipelines uses a sparse Zarr layer today (dense + `X`/layers only; the sparse `rare`-stratum genotype is read from PGEN + directly, never persisted to Zarr). Revisit if/when a real sparse Zarr + layer is added. ### Changed diff --git a/README.md b/README.md index 511ccde..aec2216 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,11 @@ + +

cellink

+ + [![Build](https://github.com/theislab/cellink/actions/workflows/build.yaml/badge.svg)](https://github.com/theislab/cellink/actions/workflows/build.yaml/badge.svg) [![License](https://img.shields.io/github/license/theislab/cellink)](https://opensource.org/licenses/Apache2.0) [![Read the Docs](https://img.shields.io/readthedocs/cellink/latest.svg?label=Read%20the%20Docs)](https://cellink-docs.readthedocs.io/) diff --git a/docs/api/tl.md b/docs/api/tl.md index 93072f4..1129d84 100644 --- a/docs/api/tl.md +++ b/docs/api/tl.md @@ -19,4 +19,6 @@ tl.beta_weighting tl.subset_genomic_region tl.subset_gene + tl.coloc_abf + tl.coloc_susie ``` diff --git a/docs/api/tl_external.md b/docs/api/tl_external.md index 819a701..7f871fd 100644 --- a/docs/api/tl_external.md +++ b/docs/api/tl_external.md @@ -2,7 +2,7 @@ Wrappers around external genetics tools. Most of these shell out to a command-line tool that has to be installed separately (e.g. PLINK, MAGMA, LDSC, TensorQTL, -SAIGE-QTL) — see each tool's own documentation for installation instructions. +SAIGE-QTL). See each tool's own documentation for installation instructions. ### eQTL mapping (jaxQTL / TensorQTL) diff --git a/docs/contributing.md b/docs/contributing.md index a170667..1f45302 100644 --- a/docs/contributing.md +++ b/docs/contributing.md @@ -52,7 +52,7 @@ When you're touching one of these wrappers, install just the extra(s) you need, `pip install -e ".[pgen,test]"`. A few (MAGMA, LDSC, PLINK/PLINK2, TensorQTL's CLI mode, SAIGE-QTL) additionally shell -out to an external binary that isn't distributed on PyPI/conda at all — tests exercising +out to an external binary that isn't distributed on PyPI/conda at all. Tests exercising those either skip via `pytest.importorskip`/an explicit binary check, or aren't run in CI and are only expected to pass in an environment where the tool is installed separately. Don't assume the base test env in CI covers every `tl.external` module. @@ -96,7 +96,7 @@ Consider enabling this option for [ruff][ruff-editors] and [prettier][prettier-e ## Writing tests This package uses [pytest][] for automated testing. Please add a test for every -function you add or change — see scanpy's {doc}`scanpy:dev/testing` guide for general +function you add or change. See scanpy's {doc}`scanpy:dev/testing` guide for general pytest conventions. Small synthetic-data tests (`tests/conftest.py` has `adata`/`gdata` fixtures built from `cellink._core.dummy_data`) are preferred over ones that depend on `tests/data/*` fixtures or network access; add a real assertion, not just a call that @@ -197,18 +197,18 @@ Public functions/classes are exported from each submodule's `__init__.py` (and l its `__all__`), and separately listed again in the matching `docs/api/*.md` page (`donordata.md`, `pp.md`, `io.md`, `tl.md`, `tl_external.md`, `pl.md`, `ml.md`, `at.md`, `utils.md`, `resources.md`, `cli.md`). These are two separate, manually-kept-in-sync -lists — nothing enforces that they match. When you add, rename, or remove a public +lists. Nothing enforces that they match. When you add, rename, or remove a public function: 1. Add it to the submodule's `__init__.py` import and `__all__` (or remove it from both). 2. Add/update the corresponding `docs/api/*.md` entry. -3. If it's a new tutorial notebook, add it to `docs/tutorials/index.md`'s toctree — a +3. If it's a new tutorial notebook, add it to `docs/tutorials/index.md`'s toctree. A notebook that exists under `docs/tutorials/` but isn't listed there won't show up in the rendered docs navigation. `ruff`'s `RUF100`/`F822` checks catch a stale `__all__` (a name listed that no longer exists), but nothing currently catches the reverse (a public name that exists but was -never added to `__all__`/the API docs) or a tutorial notebook missing from the toctree — +never added to `__all__`/the API docs) or a tutorial notebook missing from the toctree, worth double-checking by eye. [sphinx]: https://www.sphinx-doc.org/en/master/ @@ -222,7 +222,7 @@ worth double-checking by eye. The documentation is set-up to render jupyter notebooks stored in the `docs/tutorials` directory using [myst-nb][]. Currently, only notebooks in `.ipynb` format are supported that will be included with both their input and output cells. -It is your responsibility to update and re-run the notebook whenever necessary — notebooks aren't executed as +It is your responsibility to update and re-run the notebook whenever necessary. Notebooks aren't executed as part of CI, so a stale output cell showing an old error/result won't be caught automatically. If you are interested in automatically running notebooks as part of the continuous integration, diff --git a/docs/index.md b/docs/index.md index c5ce7ee..a6c1b69 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,4 +1,15 @@ ```{include} ../README.md +:end-before: '' +``` + +```{image} _static/img/cellink-logo.svg +:align: center +:width: 360px +:alt: cellink +``` + +```{include} ../README.md +:start-after: '' :end-before: '' ``` diff --git a/docs/tutorials/burden_testing.ipynb b/docs/tutorials/burden_testing.ipynb index 43eb4d3..170788f 100644 --- a/docs/tutorials/burden_testing.ipynb +++ b/docs/tutorials/burden_testing.ipynb @@ -8721,7 +8721,7 @@ "source": [ "## Step 2: Burden Testing\n", "\n", - "We now perform burden tests, which evaluate whether the **aggregate effect** of rare variants within a gene is associated with a phenotype — in this case, gene expression in CD8 Naive cells. We use different weighting schemes allow us to prioritize variants differently based on functional scores or allele frequency, including:\n", + "We now perform burden tests, which evaluate whether the **aggregate effect** of rare variants within a gene is associated with a phenotype, in this case gene expression in CD8 Naive cells. We use different weighting schemes allow us to prioritize variants differently based on functional scores or allele frequency, including:\n", "- **CADD_RAW**: Raw CADD scores.\n", "- **maf_beta**: Beta weighting based on minor allele frequency (MAF).\n", "- **tss_distance**: Distance to the transcription start site (TSS).\n", @@ -8746,7 +8746,7 @@ "metadata": {}, "source": [ "### Filtering for Rare Variants\n", - "Burden tests focus on rare variants—typically those with minor allele frequency (MAF) below 0.01 or 0.05. We filter our genotype data accordingly to isolate variants in this frequency range." + "Burden tests focus on rare variants, typically those with minor allele frequency (MAF) below 0.01 or 0.05. We filter our genotype data accordingly to isolate variants in this frequency range." ] }, { diff --git a/docs/tutorials/cell_level_ldsc_analysis.ipynb b/docs/tutorials/cell_level_ldsc_analysis.ipynb index def6186..502fa0c 100644 --- a/docs/tutorials/cell_level_ldsc_analysis.ipynb +++ b/docs/tutorials/cell_level_ldsc_analysis.ipynb @@ -207,7 +207,7 @@ ], "source": [ "# Load the dummy OneK1K dataset (genotype + single-cell expression).\n", - "# Downloaded and cached automatically — see get_dummy_onek1k docs.\n", + "# Downloaded and cached automatically. See get_dummy_onek1k docs.\n", "dd = get_dummy_onek1k()\n", "print(f\"Dataset shape: {dd.shape}\")" ] @@ -970,7 +970,7 @@ "import cellink\n", "from pathlib import Path as _Path\n", "\n", - "# Configure for your environment — edit the YAML or pass config_dict directly.\n", + "# Configure for your environment: edit the YAML or pass config_dict directly.\n", "# Configs for Docker / Singularity / local are in cellink/tl/external/config/.\n", "_ldsc_config = _Path(cellink.__file__).parent / \"tl\" / \"external\" / \"config\" / \"ldsc_singularity.yaml\"\n", "runner = configure_ldsc_runner(config_path=str(_ldsc_config))" @@ -981,7 +981,7 @@ "metadata": {}, "source": [ "#### Step 3: Prepare Data for Analysis\n", - "To speed up computation in this tutorial, we filter the data to a random subset of SNPs per chromosome. Note: in a real analysis, you would process all chromosomes without subsetting. (We still process every chromosome here — `ldsc.py --h2-cts` requires LD scores for the whole genome — just with far fewer SNPs per chromosome than a real analysis.)\n" + "To speed up computation in this tutorial, we filter the data to a random subset of SNPs per chromosome. Note: in a real analysis, you would process all chromosomes without subsetting. (We still process every chromosome here, since `ldsc.py --h2-cts` requires LD scores for the whole genome, just with far fewer SNPs per chromosome than a real analysis.)\n" ] }, { @@ -1112,14 +1112,14 @@ "\n", "We create annotation files that mark which SNPs are near cell-type-specific genes. `cellink` provides two annotation types that produce identical output formats and are fully interchangeable downstream:\n", "\n", - "- **Binary** (`make_annot_from_donor_data` / `make_annot_from_bimfile`) — SNPs within the window of a top-10% gene get `ANNOT=1`, all others `0`. Requires a `.GeneSet` file.\n", - "- **Continuous** (`make_annot_from_donor_data` / `make_annot_from_bimfile`) — each SNP gets the specificity score of the highest-scoring gene whose window overlaps it (or `0` if none). Preserves the full quantitative signal.\n", + "- **Binary** (`make_annot_from_donor_data` / `make_annot_from_bimfile`): SNPs within the window of a top-10% gene get `ANNOT=1`, all others `0`. Requires a `.GeneSet` file.\n", + "- **Continuous** (`make_annot_from_donor_data` / `make_annot_from_bimfile`): each SNP gets the specificity score of the highest-scoring gene whose window overlaps it (or `0` if none). Preserves the full quantitative signal.\n", "\n", - "**When to use which:** Binary annotations are the standard, widely-used approach (e.g. Finucane et al. LDSC-SEG) — simple to interpret, easy to compare across studies, and the right default when you just want to test whether heritability is enriched in the top N% most cell-type-specific genes. Continuous annotations avoid the arbitrary top-N% cutoff and can be more powerful when specificity varies smoothly across genes rather than splitting cleanly into high/low — they are most useful when you already have a continuous score (e.g. from SEISMIC or Duncan et al.) that you don't want to threshold away. **This tutorial uses binary annotations as the main, recommended workflow**; the continuous variant is shown afterward as an alternative for when a continuous score is preferable.\n", + "**When to use which:** Binary annotations are the standard, widely-used approach (e.g. Finucane et al. LDSC-SEG), simple to interpret, easy to compare across studies, and the right default when you just want to test whether heritability is enriched in the top N% most cell-type-specific genes. Continuous annotations avoid the arbitrary top-N% cutoff and can be more powerful when specificity varies smoothly across genes rather than splitting cleanly into high/low; they are most useful when you already have a continuous score (e.g. from SEISMIC or Duncan et al.) that you don't want to threshold away. **This tutorial uses binary annotations as the main, recommended workflow**; the continuous variant is shown afterward as an alternative for when a continuous score is preferable.\n", "\n", "Both write `CHR / BP / SNP / CM / ANNOT` files, so the downstream LD-score call is identical for both.\n", "\n", - "**Binary annotations from your own genotype data** (top-10% specificity genes → 1), using `make_annot_from_donor_data` — illustrative only; we don't carry this particular version forward (see why below):" + "**Binary annotations from your own genotype data** (top-10% specificity genes → 1), using `make_annot_from_donor_data`, illustrative only; we don't carry this particular version forward (see why below):" ] }, { @@ -1155,7 +1155,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "**Why we switch to a reference panel for the rest of this tutorial:** OneK1K's own genotype SNP IDs are positional (e.g. `1_916119_A_G`), but real GWAS summary statistics are keyed by rsID. LD scores computed from `dd`'s own genotypes can never be merged with a real GWAS downstream — and this is also how real S-LDSC analyses are actually run: LD is computed from a **reference panel** (e.g. 1000 Genomes), while the study's own single-cell data is only used to derive the specificity scores that define the annotation. So from here on we create annotations against 1000 Genomes plink files instead of `dd`'s genotypes — same functions, just swap `make_annot_from_donor_data` for `make_annot_from_bimfile`:" + "**Why we switch to a reference panel for the rest of this tutorial:** OneK1K's own genotype SNP IDs are positional (e.g. `1_916119_A_G`), but real GWAS summary statistics are keyed by rsID. LD scores computed from `dd`'s own genotypes can never be merged with a real GWAS downstream. This is also how real S-LDSC analyses are actually run: LD is computed from a **reference panel** (e.g. 1000 Genomes), while the study's own single-cell data is only used to derive the specificity scores that define the annotation. So from here on we create annotations against 1000 Genomes plink files instead of `dd`'s genotypes, same functions, just swap `make_annot_from_donor_data` for `make_annot_from_bimfile`:" ] }, { @@ -1212,7 +1212,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "**Continuous annotations (alternative)** — instead of thresholding to the top 10%, assign each SNP the actual specificity score of its nearest gene window. Use this if you'd rather not pick an arbitrary cutoff. Shown first from your own donor data (illustrative, same reasoning as above), then against the 1000G reference panel (the version this tutorial actually uses downstream):" + "**Continuous annotations (alternative)**: instead of thresholding to the top 10%, assign each SNP the actual specificity score of its nearest gene window. Use this if you'd rather not pick an arbitrary cutoff. Shown first from your own donor data (illustrative, same reasoning as above), then against the 1000G reference panel (the version this tutorial actually uses downstream):" ] }, { @@ -1252,7 +1252,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Against the 1000 Genomes reference panel (no `DonorData` needed) — useful when running annotations independently from genotype data, and the version that's rsID-consistent with real GWAS:" + "Against the 1000 Genomes reference panel (no `DonorData` needed), useful when running annotations independently from genotype data, and the version that's rsID-consistent with real GWAS:" ] }, { @@ -1286,7 +1286,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "**Comparing binary vs continuous** — both files have identical columns; the only difference is the ANNOT distribution. Comparing the 1000G-based versions, since those are the ones used downstream:" + "**Comparing binary vs continuous**: both files have identical columns; the only difference is the ANNOT distribution. Comparing the 1000G-based versions, since those are the ones used downstream:" ] }, { @@ -1384,7 +1384,7 @@ "metadata": {}, "source": [ "#### Step 6: Prepare Reference LD Scores and Weights\n", - "For the final analysis, we need baseline LD scores and regression weights. These control for genomic confounders and ensure proper statistical inference. We download these from the 1000 Genomes reference panel — the same real, rsID-based files used in published S-LDSC analyses, so they merge correctly with the real GWAS downloaded in Step 7." + "For the final analysis, we need baseline LD scores and regression weights. These control for genomic confounders and ensure proper statistical inference. We download these from the 1000 Genomes reference panel, the same real, rsID-based files used in published S-LDSC analyses, so they merge correctly with the real GWAS downloaded in Step 7." ] }, { @@ -1609,7 +1609,7 @@ "metadata": {}, "source": [ "#### Step 8: Run Cell-Type-Specific Heritability Analysis\n", - "Finally, we run the cell-type-specific heritability analysis, which tests whether each cell type shows enrichment for trait heritability. `ref_ld_chr`/`w_ld_chr` use the real 1000G baseline/weights downloaded in Step 6 — the same rsID convention as the GWAS sumstats and the cell-type-specific LD scores, so the merge step finds real overlap." + "Finally, we run the cell-type-specific heritability analysis, which tests whether each cell type shows enrichment for trait heritability. `ref_ld_chr`/`w_ld_chr` use the real 1000G baseline/weights downloaded in Step 6, the same rsID convention as the GWAS sumstats and the cell-type-specific LD scores, so the merge step finds real overlap." ] }, { diff --git a/docs/tutorials/colocalization.ipynb b/docs/tutorials/colocalization.ipynb index 73164dd..1e54c11 100644 --- a/docs/tutorials/colocalization.ipynb +++ b/docs/tutorials/colocalization.ipynb @@ -1139,6 +1139,97 @@ " print(\"Inconclusive\")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Native Python Implementation (No R Dependency)\n", + "\n", + "`cellink.tl.coloc_abf` implements the same Approximate Bayes Factor statistic (Giambartolomei et al. 2014) directly in NumPy, thus no R installation or `rpy2` is required. It is a drop-in numerical equivalent of `coloc::coloc.abf()`. We verify this below against the real R result computed above, on this tutorial's own data, and separately against 6 independent synthetic test cases (max absolute difference 8.3e-15 across every PP0-PP4 value, which is a floating-point precision, not a real numerical discrepancy)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from cellink.tl import DEFAULT_PRIOR_VAR_CC, coloc_abf\n", + "\n", + "native_pp = coloc_abf(\n", + " beta1=eqtl_harm[\"beta\"].to_numpy(),\n", + " se1=np.sqrt(eqtl_harm[\"varbeta\"].to_numpy()),\n", + " beta2=gwas_harm[\"beta\"].to_numpy(),\n", + " se2=np.sqrt(gwas_harm[\"varbeta\"].to_numpy()),\n", + " p1=coloc_p1,\n", + " p2=coloc_p2,\n", + " p12=coloc_p12,\n", + " prior_var2=DEFAULT_PRIOR_VAR_CC,\n", + ")\n", + "print(\"Native Python coloc_abf:\")\n", + "for k in [\"PP0\", \"PP1\", \"PP2\", \"PP3\", \"PP4\"]:\n", + " print(f\" {k}: {native_pp[k]:.6f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Does the native result match R exactly, on this exact data?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "r_pp = {\n", + " \"PP0\": coloc_result[\"PP.H0.abf\"],\n", + " \"PP1\": coloc_result[\"PP.H1.abf\"],\n", + " \"PP2\": coloc_result[\"PP.H2.abf\"],\n", + " \"PP3\": coloc_result[\"PP.H3.abf\"],\n", + " \"PP4\": coloc_result[\"PP.H4.abf\"],\n", + "}\n", + "print(f\"{'':6}{'R (rpy2)':>12}{'Native':>12}{'abs diff':>14}\")\n", + "max_diff = 0.0\n", + "for k in r_pp:\n", + " diff = abs(r_pp[k] - native_pp[k])\n", + " max_diff = max(max_diff, diff)\n", + " print(f\"{k:6}{r_pp[k]:12.6f}{native_pp[k]:12.6f}{diff:14.2e}\")\n", + "print(f\"\\nMax absolute difference: {max_diff:.2e}\")\n", + "assert max_diff < 1e-6, \"native coloc_abf disagrees with R coloc::coloc.abf() beyond floating-point noise\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### `coloc_susie` in practice: a self-contained demonstration\n", + "\n", + "Uses synthetic genotypes/phenotypes (not this tutorial's real eQTL/GWAS data, which only has one credible set per side and so doesn't exercise the multi-signal case `coloc_susie` is for) to show the actual calling convention: fit SuSiE on both sides, extract `lbf_variable` and each fit's own credible-set indices, then colocalize." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from cellink.tl import coloc_susie\n", + "\n", + "n_snps = 30\n", + "lbf1 = np.zeros((2, n_snps))\n", + "lbf1[0, 5] = 25.0\n", + "lbf1[1, 20] = 22.0\n", + "lbf2 = np.zeros((1, n_snps))\n", + "lbf2[0, 5] = 18.0\n", + "\n", + "res = coloc_susie(lbf1, lbf2, cs1_index=[0, 1], cs2_index=[0])\n", + "print(res[[\"idx1\", \"idx2\", \"hit1\", \"hit2\", \"PP0\", \"PP1\", \"PP2\", \"PP3\", \"PP4\"]])" + ] + }, { "cell_type": "markdown", "metadata": {}, diff --git a/docs/tutorials/donordata_basics.ipynb b/docs/tutorials/donordata_basics.ipynb index 07e4167..273170d 100644 --- a/docs/tutorials/donordata_basics.ipynb +++ b/docs/tutorials/donordata_basics.ipynb @@ -62,7 +62,7 @@ { "data": { "text/plain": [ - "AnnData object with n_obs × n_vars = 6 × 4\n", + "AnnData object with n_obs \u00d7 n_vars = 6 \u00d7 4\n", " obs: 'donor_id', 'celltype'" ] }, @@ -230,25 +230,25 @@ { "data": { "text/html": [ - "
╔═ DonorData(n_donors=3, n_cells_per_donor=[2-2], donor_id='donor_id') ═════════════════════════╗\n",
-       "║ ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ║\n",
-       "║ ┃ G (donors)                                  ┃ C (cells)                                   ┃ ║\n",
-       "║ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ ║\n",
-       "║ │ AnnData object with n_obs × n_vars = 3 × 5  │ AnnData object with n_obs × n_vars = 6 × 4  │ ║\n",
-       "║ │     var: 'chrom', 'pos'                     │     obs: 'donor_id', 'celltype'             │ ║\n",
-       "║ └─────────────────────────────────────────────┴─────────────────────────────────────────────┘ ║\n",
-       "╚═══════════════════════════════════════════════════════════════════════════════════════════════╝\n",
+       "
\u2554\u2550 DonorData(n_donors=3, n_cells_per_donor=[2-2], donor_id='donor_id') \u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2557\n",
+       "\u2551 \u250f\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2533\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2513 \u2551\n",
+       "\u2551 \u2503 G (donors)                                  \u2503 C (cells)                                   \u2503 \u2551\n",
+       "\u2551 \u2521\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2547\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2529 \u2551\n",
+       "\u2551 \u2502 AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 5  \u2502 AnnData object with n_obs \u00d7 n_vars = 6 \u00d7 4  \u2502 \u2551\n",
+       "\u2551 \u2502     var: 'chrom', 'pos'                     \u2502     obs: 'donor_id', 'celltype'             \u2502 \u2551\n",
+       "\u2551 \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2551\n",
+       "\u255a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u255d\n",
        "
\n"
       ],
       "text/plain": [
-       "╔═\u001b[1;38;5;197m DonorData(n_donors=3, n_cells_per_donor=[2-2], donor_id='donor_id') \u001b[0m═════════════════════════╗\n",
-       "║ ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ║\n",
-       "║ ┃\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mG (donors)                                 \u001b[0m\u001b[1;38;5;197m \u001b[0m┃\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mC (cells)                                  \u001b[0m\u001b[1;38;5;197m \u001b[0m┃ ║\n",
-       "║ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ ║\n",
-       "║ │ AnnData object with n_obs × n_vars = 3 × 5  │ AnnData object with n_obs × n_vars = 6 × 4  │ ║\n",
-       "║ │     var: 'chrom', 'pos'                     │     obs: \u001b[1;38;5;197m'donor_id', \u001b[0m'celltype'             │ ║\n",
-       "║ └─────────────────────────────────────────────┴─────────────────────────────────────────────┘ ║\n",
-       "╚═══════════════════════════════════════════════════════════════════════════════════════════════╝\n"
+       "\u2554\u2550\u001b[1;38;5;197m DonorData(n_donors=3, n_cells_per_donor=[2-2], donor_id='donor_id') \u001b[0m\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2557\n",
+       "\u2551 \u250f\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2533\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2513 \u2551\n",
+       "\u2551 \u2503\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mG (donors)                                 \u001b[0m\u001b[1;38;5;197m \u001b[0m\u2503\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mC (cells)                                  \u001b[0m\u001b[1;38;5;197m \u001b[0m\u2503 \u2551\n",
+       "\u2551 \u2521\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2547\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2529 \u2551\n",
+       "\u2551 \u2502 AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 5  \u2502 AnnData object with n_obs \u00d7 n_vars = 6 \u00d7 4  \u2502 \u2551\n",
+       "\u2551 \u2502     var: 'chrom', 'pos'                     \u2502     obs: \u001b[1;38;5;197m'donor_id', \u001b[0m'celltype'             \u2502 \u2551\n",
+       "\u2551 \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2551\n",
+       "\u255a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u255d\n"
       ]
      },
      "metadata": {},
@@ -423,7 +423,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[2026-07-11 17:36:22,567] INFO:cellink.io._pgen: Total: 100 samples × 1,000 variants (dense)\n"
+      "[2026-07-11 17:36:22,567] INFO:cellink.io._pgen: Total: 100 samples \u00d7 1,000 variants (dense)\n"
      ]
     },
     {
@@ -437,7 +437,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[2026-07-11 17:36:22,607] INFO:cellink.io._pgen: Creating dense X (100 × 1,000) chunks=(50, 200) ...\n"
+      "[2026-07-11 17:36:22,607] INFO:cellink.io._pgen: Creating dense X (100 \u00d7 1,000) chunks=(50, 200) ...\n"
      ]
     },
     {
@@ -460,7 +460,7 @@
      "output_type": "stream",
      "text": [
       "\r",
-      "  simulated_genotype_calls:  40%|████      | 400/1000 [00:00<00:00, 2444.58var/s]"
+      "  simulated_genotype_calls:  40%|\u2588\u2588\u2588\u2588      | 400/1000 [00:00<00:00, 2444.58var/s]"
      ]
     },
     {
@@ -468,7 +468,7 @@
      "output_type": "stream",
      "text": [
       "\r",
-      "  simulated_genotype_calls: 100%|██████████| 1000/1000 [00:00<00:00, 5834.68var/s]"
+      "  simulated_genotype_calls: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 1000/1000 [00:00<00:00, 5834.68var/s]"
      ]
     },
     {
@@ -489,7 +489,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[2026-07-11 17:36:22,968] INFO:cellink.io._pgen: ✓ Done → /tmp/tmpwu4s7bxq/pgen_demo.zarr\n"
+      "[2026-07-11 17:36:22,968] INFO:cellink.io._pgen: \u2713 Done \u2192 /tmp/tmpwu4s7bxq/pgen_demo.zarr\n"
      ]
     },
     {
@@ -684,25 +684,25 @@
     {
      "data": {
       "text/html": [
-       "
╔═ DonorData(n_donors=2, n_cells_per_donor=[2-2], donor_id='donor_id') ═══════════════════════════════════════╗\n",
-       "║ ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ║\n",
-       "║ ┃ G (donors)                                         ┃ C (cells)                                          ┃ ║\n",
-       "║ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ ║\n",
-       "║ │ View of AnnData object with n_obs × n_vars = 2 × 5 │ View of AnnData object with n_obs × n_vars = 4 × 4 │ ║\n",
-       "║ │     var: 'chrom', 'pos'                            │     obs: 'donor_id', 'celltype'                    │ ║\n",
-       "║ └────────────────────────────────────────────────────┴────────────────────────────────────────────────────┘ ║\n",
-       "╚═════════════════════════════════════════════════════════════════════════════════════════════════════════════╝\n",
+       "
\u2554\u2550 DonorData(n_donors=2, n_cells_per_donor=[2-2], donor_id='donor_id') \u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2557\n",
+       "\u2551 \u250f\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2533\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2513 \u2551\n",
+       "\u2551 \u2503 G (donors)                                         \u2503 C (cells)                                          \u2503 \u2551\n",
+       "\u2551 \u2521\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2547\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2529 \u2551\n",
+       "\u2551 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 2 \u00d7 5 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 4 \u00d7 4 \u2502 \u2551\n",
+       "\u2551 \u2502     var: 'chrom', 'pos'                            \u2502     obs: 'donor_id', 'celltype'                    \u2502 \u2551\n",
+       "\u2551 \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2551\n",
+       "\u255a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u255d\n",
        "
\n"
       ],
       "text/plain": [
-       "╔═\u001b[1;38;5;197m DonorData(n_donors=2, n_cells_per_donor=[2-2], donor_id='donor_id') \u001b[0m═══════════════════════════════════════╗\n",
-       "║ ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ║\n",
-       "║ ┃\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mG (donors)                                        \u001b[0m\u001b[1;38;5;197m \u001b[0m┃\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mC (cells)                                         \u001b[0m\u001b[1;38;5;197m \u001b[0m┃ ║\n",
-       "║ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ ║\n",
-       "║ │ View of AnnData object with n_obs × n_vars = 2 × 5 │ View of AnnData object with n_obs × n_vars = 4 × 4 │ ║\n",
-       "║ │     var: 'chrom', 'pos'                            │     obs: \u001b[1;38;5;197m'donor_id', \u001b[0m'celltype'                    │ ║\n",
-       "║ └────────────────────────────────────────────────────┴────────────────────────────────────────────────────┘ ║\n",
-       "╚═════════════════════════════════════════════════════════════════════════════════════════════════════════════╝\n"
+       "\u2554\u2550\u001b[1;38;5;197m DonorData(n_donors=2, n_cells_per_donor=[2-2], donor_id='donor_id') \u001b[0m\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2557\n",
+       "\u2551 \u250f\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2533\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2513 \u2551\n",
+       "\u2551 \u2503\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mG (donors)                                        \u001b[0m\u001b[1;38;5;197m \u001b[0m\u2503\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mC (cells)                                         \u001b[0m\u001b[1;38;5;197m \u001b[0m\u2503 \u2551\n",
+       "\u2551 \u2521\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2547\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2529 \u2551\n",
+       "\u2551 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 2 \u00d7 5 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 4 \u00d7 4 \u2502 \u2551\n",
+       "\u2551 \u2502     var: 'chrom', 'pos'                            \u2502     obs: \u001b[1;38;5;197m'donor_id', \u001b[0m'celltype'                    \u2502 \u2551\n",
+       "\u2551 \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2551\n",
+       "\u255a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u255d\n"
       ]
      },
      "metadata": {},
@@ -737,25 +737,25 @@
     {
      "data": {
       "text/html": [
-       "
╔═ DonorData(n_donors=3, n_cells_per_donor=[1-1], donor_id='donor_id') ═══════════════════════════════════════╗\n",
-       "║ ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ║\n",
-       "║ ┃ G (donors)                                         ┃ C (cells)                                          ┃ ║\n",
-       "║ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ ║\n",
-       "║ │ View of AnnData object with n_obs × n_vars = 3 × 5 │ View of AnnData object with n_obs × n_vars = 3 × 4 │ ║\n",
-       "║ │     var: 'chrom', 'pos'                            │     obs: 'donor_id', 'celltype'                    │ ║\n",
-       "║ └────────────────────────────────────────────────────┴────────────────────────────────────────────────────┘ ║\n",
-       "╚═════════════════════════════════════════════════════════════════════════════════════════════════════════════╝\n",
+       "
\u2554\u2550 DonorData(n_donors=3, n_cells_per_donor=[1-1], donor_id='donor_id') \u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2557\n",
+       "\u2551 \u250f\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2533\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2513 \u2551\n",
+       "\u2551 \u2503 G (donors)                                         \u2503 C (cells)                                          \u2503 \u2551\n",
+       "\u2551 \u2521\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2547\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2529 \u2551\n",
+       "\u2551 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 5 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 4 \u2502 \u2551\n",
+       "\u2551 \u2502     var: 'chrom', 'pos'                            \u2502     obs: 'donor_id', 'celltype'                    \u2502 \u2551\n",
+       "\u2551 \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2551\n",
+       "\u255a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u255d\n",
        "
\n"
       ],
       "text/plain": [
-       "╔═\u001b[1;38;5;197m DonorData(n_donors=3, n_cells_per_donor=[1-1], donor_id='donor_id') \u001b[0m═══════════════════════════════════════╗\n",
-       "║ ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ║\n",
-       "║ ┃\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mG (donors)                                        \u001b[0m\u001b[1;38;5;197m \u001b[0m┃\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mC (cells)                                         \u001b[0m\u001b[1;38;5;197m \u001b[0m┃ ║\n",
-       "║ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ ║\n",
-       "║ │ View of AnnData object with n_obs × n_vars = 3 × 5 │ View of AnnData object with n_obs × n_vars = 3 × 4 │ ║\n",
-       "║ │     var: 'chrom', 'pos'                            │     obs: \u001b[1;38;5;197m'donor_id', \u001b[0m'celltype'                    │ ║\n",
-       "║ └────────────────────────────────────────────────────┴────────────────────────────────────────────────────┘ ║\n",
-       "╚═════════════════════════════════════════════════════════════════════════════════════════════════════════════╝\n"
+       "\u2554\u2550\u001b[1;38;5;197m DonorData(n_donors=3, n_cells_per_donor=[1-1], donor_id='donor_id') \u001b[0m\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2557\n",
+       "\u2551 \u250f\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2533\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2513 \u2551\n",
+       "\u2551 \u2503\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mG (donors)                                        \u001b[0m\u001b[1;38;5;197m \u001b[0m\u2503\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mC (cells)                                         \u001b[0m\u001b[1;38;5;197m \u001b[0m\u2503 \u2551\n",
+       "\u2551 \u2521\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2547\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2529 \u2551\n",
+       "\u2551 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 5 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 4 \u2502 \u2551\n",
+       "\u2551 \u2502     var: 'chrom', 'pos'                            \u2502     obs: \u001b[1;38;5;197m'donor_id', \u001b[0m'celltype'                    \u2502 \u2551\n",
+       "\u2551 \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2551\n",
+       "\u255a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u255d\n"
       ]
      },
      "metadata": {},
@@ -799,25 +799,25 @@
     {
      "data": {
       "text/html": [
-       "
╔═ DonorData(n_donors=3, n_cells_per_donor=[1-1], donor_id='donor_id') ═══════════════════════════════════════╗\n",
-       "║ ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ║\n",
-       "║ ┃ G (donors)                                         ┃ C (cells)                                          ┃ ║\n",
-       "║ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ ║\n",
-       "║ │ View of AnnData object with n_obs × n_vars = 3 × 5 │ View of AnnData object with n_obs × n_vars = 3 × 4 │ ║\n",
-       "║ │     var: 'chrom', 'pos'                            │     obs: 'donor_id', 'celltype'                    │ ║\n",
-       "║ └────────────────────────────────────────────────────┴────────────────────────────────────────────────────┘ ║\n",
-       "╚═════════════════════════════════════════════════════════════════════════════════════════════════════════════╝\n",
+       "
\u2554\u2550 DonorData(n_donors=3, n_cells_per_donor=[1-1], donor_id='donor_id') \u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2557\n",
+       "\u2551 \u250f\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2533\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2513 \u2551\n",
+       "\u2551 \u2503 G (donors)                                         \u2503 C (cells)                                          \u2503 \u2551\n",
+       "\u2551 \u2521\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2547\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2529 \u2551\n",
+       "\u2551 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 5 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 4 \u2502 \u2551\n",
+       "\u2551 \u2502     var: 'chrom', 'pos'                            \u2502     obs: 'donor_id', 'celltype'                    \u2502 \u2551\n",
+       "\u2551 \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2551\n",
+       "\u255a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u255d\n",
        "
\n"
       ],
       "text/plain": [
-       "╔═\u001b[1;38;5;197m DonorData(n_donors=3, n_cells_per_donor=[1-1], donor_id='donor_id') \u001b[0m═══════════════════════════════════════╗\n",
-       "║ ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ║\n",
-       "║ ┃\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mG (donors)                                        \u001b[0m\u001b[1;38;5;197m \u001b[0m┃\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mC (cells)                                         \u001b[0m\u001b[1;38;5;197m \u001b[0m┃ ║\n",
-       "║ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ ║\n",
-       "║ │ View of AnnData object with n_obs × n_vars = 3 × 5 │ View of AnnData object with n_obs × n_vars = 3 × 4 │ ║\n",
-       "║ │     var: 'chrom', 'pos'                            │     obs: \u001b[1;38;5;197m'donor_id', \u001b[0m'celltype'                    │ ║\n",
-       "║ └────────────────────────────────────────────────────┴────────────────────────────────────────────────────┘ ║\n",
-       "╚═════════════════════════════════════════════════════════════════════════════════════════════════════════════╝\n"
+       "\u2554\u2550\u001b[1;38;5;197m DonorData(n_donors=3, n_cells_per_donor=[1-1], donor_id='donor_id') \u001b[0m\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2557\n",
+       "\u2551 \u250f\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2533\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2513 \u2551\n",
+       "\u2551 \u2503\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mG (donors)                                        \u001b[0m\u001b[1;38;5;197m \u001b[0m\u2503\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mC (cells)                                         \u001b[0m\u001b[1;38;5;197m \u001b[0m\u2503 \u2551\n",
+       "\u2551 \u2521\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2547\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2529 \u2551\n",
+       "\u2551 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 5 \u2502 View of AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 4 \u2502 \u2551\n",
+       "\u2551 \u2502     var: 'chrom', 'pos'                            \u2502     obs: \u001b[1;38;5;197m'donor_id', \u001b[0m'celltype'                    \u2502 \u2551\n",
+       "\u2551 \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2551\n",
+       "\u255a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u255d\n"
       ]
      },
      "metadata": {},
@@ -965,27 +965,27 @@
     {
      "data": {
       "text/html": [
-       "
╔═ DonorData(n_donors=3, n_cells_per_donor=[2-2], donor_id='donor_id') ═════════════════════════╗\n",
-       "║ ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ║\n",
-       "║ ┃ G (donors)                                  ┃ C (cells)                                   ┃ ║\n",
-       "║ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ ║\n",
-       "║ │ AnnData object with n_obs × n_vars = 3 × 5  │ AnnData object with n_obs × n_vars = 6 × 4  │ ║\n",
-       "║ │     var: 'chrom', 'pos'                     │     obs: 'donor_id', 'celltype'             │ ║\n",
-       "║ │     obsm: 'mean_expr'                       │                                             │ ║\n",
-       "║ └─────────────────────────────────────────────┴─────────────────────────────────────────────┘ ║\n",
-       "╚═══════════════════════════════════════════════════════════════════════════════════════════════╝\n",
+       "
\u2554\u2550 DonorData(n_donors=3, n_cells_per_donor=[2-2], donor_id='donor_id') \u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2557\n",
+       "\u2551 \u250f\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2533\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2513 \u2551\n",
+       "\u2551 \u2503 G (donors)                                  \u2503 C (cells)                                   \u2503 \u2551\n",
+       "\u2551 \u2521\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2547\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2529 \u2551\n",
+       "\u2551 \u2502 AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 5  \u2502 AnnData object with n_obs \u00d7 n_vars = 6 \u00d7 4  \u2502 \u2551\n",
+       "\u2551 \u2502     var: 'chrom', 'pos'                     \u2502     obs: 'donor_id', 'celltype'             \u2502 \u2551\n",
+       "\u2551 \u2502     obsm: 'mean_expr'                       \u2502                                             \u2502 \u2551\n",
+       "\u2551 \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2551\n",
+       "\u255a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u255d\n",
        "
\n"
       ],
       "text/plain": [
-       "╔═\u001b[1;38;5;197m DonorData(n_donors=3, n_cells_per_donor=[2-2], donor_id='donor_id') \u001b[0m═════════════════════════╗\n",
-       "║ ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ║\n",
-       "║ ┃\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mG (donors)                                 \u001b[0m\u001b[1;38;5;197m \u001b[0m┃\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mC (cells)                                  \u001b[0m\u001b[1;38;5;197m \u001b[0m┃ ║\n",
-       "║ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ ║\n",
-       "║ │ AnnData object with n_obs × n_vars = 3 × 5  │ AnnData object with n_obs × n_vars = 6 × 4  │ ║\n",
-       "║ │     var: 'chrom', 'pos'                     │     obs: \u001b[1;38;5;197m'donor_id', \u001b[0m'celltype'             │ ║\n",
-       "║ │     obsm: 'mean_expr'                       │                                             │ ║\n",
-       "║ └─────────────────────────────────────────────┴─────────────────────────────────────────────┘ ║\n",
-       "╚═══════════════════════════════════════════════════════════════════════════════════════════════╝\n"
+       "\u2554\u2550\u001b[1;38;5;197m DonorData(n_donors=3, n_cells_per_donor=[2-2], donor_id='donor_id') \u001b[0m\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2557\n",
+       "\u2551 \u250f\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2533\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2513 \u2551\n",
+       "\u2551 \u2503\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mG (donors)                                 \u001b[0m\u001b[1;38;5;197m \u001b[0m\u2503\u001b[1;38;5;197m \u001b[0m\u001b[1;38;5;197mC (cells)                                  \u001b[0m\u001b[1;38;5;197m \u001b[0m\u2503 \u2551\n",
+       "\u2551 \u2521\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2547\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2501\u2529 \u2551\n",
+       "\u2551 \u2502 AnnData object with n_obs \u00d7 n_vars = 3 \u00d7 5  \u2502 AnnData object with n_obs \u00d7 n_vars = 6 \u00d7 4  \u2502 \u2551\n",
+       "\u2551 \u2502     var: 'chrom', 'pos'                     \u2502     obs: \u001b[1;38;5;197m'donor_id', \u001b[0m'celltype'             \u2502 \u2551\n",
+       "\u2551 \u2502     obsm: 'mean_expr'                       \u2502                                             \u2502 \u2551\n",
+       "\u2551 \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2534\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2551\n",
+       "\u255a\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u255d\n"
       ]
      },
      "metadata": {},
@@ -1013,6 +1013,41 @@
     "dd_loaded"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "zarr-chunking-md",
+   "metadata": {},
+   "source": [
+    "`write_zarr_dd`/`read_zarr_dd` save the same `G`/`C` pair to a Zarr v3 store\n",
+    "instead of a single HDF5 file. The read side is lazy: a dense `X` (or a dense\n",
+    "layer) comes back Dask-backed, so selecting a subset (one cell type of `C`,\n",
+    "say) never touches the rest of the data on disk.\n",
+    "\n",
+    "On the write side, a dense `X`/layer gets sane chunk sizes automatically: a\n",
+    "Dask array keeps its own existing chunks (via `dask.array.to_zarr`, which\n",
+    "preserves them exactly), and a plain array gets chunks capped at 4096 per\n",
+    "axis. Pass `x_chunks` to `write_zarr_dd` to override this: a single\n",
+    "`(rows, cols)` tuple applies to both `G.X` and `C.X`, or a `{\"G\": ..., \"C\": ...}`\n",
+    "dict sets them independently. This matters at genome/transcriptome scale: an\n",
+    "auto-chosen chunk shape that doesn't respect how the data is actually laid out\n",
+    "can make both the write and every later read far slower than necessary."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "zarr-chunking-code",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from cellink.io import read_zarr_dd\n",
+    "\n",
+    "dd.write_zarr_dd(str(tmp_dir / \"dd_basics.dd.zarr\"), x_chunks=(2, 2))\n",
+    "\n",
+    "dd_loaded_zarr = read_zarr_dd(str(tmp_dir / \"dd_basics.dd.zarr\"))\n",
+    "dd_loaded_zarr.C.X  # Dask-backed; call .compute() to materialize"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "aeaa9a62",
diff --git a/docs/tutorials/ehrdataset.ipynb b/docs/tutorials/ehrdataset.ipynb
index 81ef131..aa4ae12 100644
--- a/docs/tutorials/ehrdataset.ipynb
+++ b/docs/tutorials/ehrdataset.ipynb
@@ -8,7 +8,7 @@
     "\n",
     "The `cellink` package was originally developed for integrating single-cell measurements with donor-level **genetic data**, but its modular design allows you to swap in **any structured donor-level modality**. In this tutorial, we demonstrate how to use **electronic health records (EHR)** as the donor-level input, replacing genotypes. \n",
     "\n",
-    "This unlocks a wide range of applications — for example, using blood pressure, lab values, or medical history in combination with cell-level transcriptomics."
+    "This unlocks a wide range of applications, for example using blood pressure, lab values, or medical history in combination with cell-level transcriptomics."
    ]
   },
   {
diff --git a/docs/tutorials/gsmap.ipynb b/docs/tutorials/gsmap.ipynb
index 890c65c..503c2a4 100644
--- a/docs/tutorials/gsmap.ipynb
+++ b/docs/tutorials/gsmap.ipynb
@@ -9,8 +9,8 @@
     "\n",
     "This tutorial demonstrates how to use [gsMap](https://github.com/JianYang-Lab/gsMap) to map GWAS-associated cells at single-cell spatial resolution. Two cellink helpers are used to make the execution even smoother:\n",
     "\n",
-    "- `format_gsmap_sumstats` — converts a GWAS DataFrame into the SNP/A1/A2/Z/N format gsMap expects, handling column remapping and mixed-type coercion automatically\n",
-    "- `load_gsmap_results` — reads the output CSV files back into DataFrames for downstream analysis in Python\n",
+    "- `format_gsmap_sumstats`: converts a GWAS DataFrame into the SNP/A1/A2/Z/N format gsMap expects, handling column remapping and mixed-type coercion automatically\n",
+    "- `load_gsmap_results`: reads the output CSV files back into DataFrames for downstream analysis in Python\n",
     "\n",
     "**gsMap** works by assigning a gene specificity score (GSS) to each spatial spot, linking GSS to nearby SNPs, and using stratified LD score regression (S-LDSC) to test whether spots with high gene specificity are enriched for GWAS heritability. The Cauchy combination test then aggregates spot-level p-values into region-level associations."
    ]
@@ -61,7 +61,7 @@
     "from cellink.resources import get_gwas_catalog_study_summary_stats\n",
     "from cellink.tl.external import format_gsmap_sumstats, load_gsmap_results\n",
     "\n",
-    "# Paths — adjust to your environment\n",
+    "# Paths: adjust to your environment\n",
     "GSMAP_RESOURCE_DIR = Path(\"./gsMap_resource\")\n",
     "WORKDIR = Path(\"./gsmap_output\")\n",
     "SAMPLE_NAME = \"V1_Adult_Mouse_Brain_Coronal\"\n",
@@ -181,7 +181,7 @@
     "sc.pl.spatial(\n",
     "    adata,\n",
     "    color=\"domain\",\n",
-    "    title=\"Mouse Brain Coronal — Spatial Domains (subsampled)\",\n",
+    "    title=\"Mouse Brain Coronal: Spatial Domains (subsampled)\",\n",
     ")"
    ]
   },
@@ -220,7 +220,7 @@
    "source": [
     "## 2. Prepare GWAS Summary Statistics\n",
     "\n",
-    "gsMap requires GWAS summary statistics with columns SNP, A1, A2, Z, N (gzip-compressed). We download IQ GWAS summary statistics from the GWAS Catalog (Savage et al. 2018, Nature Genetics — the same study used in the gsMap paper).\n"
+    "gsMap requires GWAS summary statistics with columns SNP, A1, A2, Z, N (gzip-compressed). We download IQ GWAS summary statistics from the GWAS Catalog (Savage et al. 2018, Nature Genetics, the same study used in the gsMap paper).\n"
    ]
   },
   {
@@ -253,7 +253,7 @@
    "source": [
     "Two things to be aware of with GWAS Catalog harmonised files:\n",
     "\n",
-    "1. **Duplicate beta columns**: this file has both `hm_beta` and a plain `beta` column containing identical values. `gsmap format_sumstats` auto-detects both as beta columns and raises `ValueError: Found 2 different BETA columns`. `format_gsmap_sumstats` avoids this by writing only the columns you explicitly name to the temporary file — so gsMap never sees the extra `beta` column.\n",
+    "1. **Duplicate beta columns**: this file has both `hm_beta` and a plain `beta` column containing identical values. `gsmap format_sumstats` auto-detects both as beta columns and raises `ValueError: Found 2 different BETA columns`. `format_gsmap_sumstats` avoids this by writing only the columns you explicitly name to the temporary file, so gsMap never sees the extra `beta` column.\n",
     "\n",
     "2. **INFO score filtering**: the file has a `mininfo` column with values around 0.6–0.8. `gsmap format_sumstats` defaults to `info_min=0.9`, which silently filters out every single SNP, leaving 0 SNPs in the output and causing the downstream jackknife regression to fail with `ValueError: More dimensions than datapoints`. Set `info_min=0.0` to skip INFO filtering for this file.\n"
    ]
@@ -331,7 +331,7 @@
    "source": [
     "## 3. Run gsMap (Quick Mode)\n",
     "\n",
-    "Quick mode runs the entire pipeline — GNN latent representations, gene specificity scores, LD score computation, spatial LDSC, and Cauchy combination — in a single command, using pre-built SNP-gene weight matrices from the resource bundle. This is the recommended starting point.\n",
+    "Quick mode runs the entire pipeline: GNN latent representations, gene specificity scores, LD score computation, spatial LDSC, and Cauchy combination, in a single command, using pre-built SNP-gene weight matrices from the resource bundle. This is the recommended starting point.\n",
     "\n",
     "The mouse Visium data uses mouse gene names, so we pass `--homolog_file` to convert them to human gene names before linking to the human GWAS.\n",
     "\n",
@@ -442,8 +442,8 @@
     "## 4. Load and Explore Results\n",
     "\n",
     "gsMap writes two key output files:\n",
-    "- **`spatial_ldsc/`** — per-spot association results (beta, se, z, p)\n",
-    "- **`cauchy_combination/`** — region-level aggregated p-values\n",
+    "- **`spatial_ldsc/`**: per-spot association results (beta, se, z, p)\n",
+    "- **`cauchy_combination/`**: region-level aggregated p-values\n",
     "\n",
     "We use `load_gsmap_results` to read both into DataFrames.\n"
    ]
@@ -687,7 +687,7 @@
     "    fig, ax = plt.subplots(figsize=(7, 5))\n",
     "    ax.barh(plot_df.index.astype(str), plot_df[\"-log10p\"], color=\"steelblue\")\n",
     "    ax.axvline(-np.log10(0.05), color=\"red\", linestyle=\"--\", label=\"p = 0.05\")\n",
-    "    ax.set_xlabel(\"-log₁₀(p) — Cauchy combination\")\n",
+    "    ax.set_xlabel(\"-log₁₀(p), Cauchy combination\")\n",
     "    ax.set_title(\"IQ-associated spatial domains (mouse brain coronal)\")\n",
     "    ax.legend()\n",
     "    plt.tight_layout()\n",
@@ -784,7 +784,7 @@
    "source": [
     "### Step-by-step mode with custom LD reference\n",
     "\n",
-    "Step-by-step mode gives full control over each pipeline stage and is useful when you want to use a custom GTF, run conditional analyses, or tune GNN parameters. It runs `generate_ldscore` for all 22 chromosomes and is compute-intensive — submit to a cluster in practice.\n",
+    "Step-by-step mode gives full control over each pipeline stage and is useful when you want to use a custom GTF, run conditional analyses, or tune GNN parameters. It runs `generate_ldscore` for all 22 chromosomes and is compute-intensive, so submit to a cluster in practice.\n",
     "\n",
     "The steps are:\n",
     "\n",
diff --git a/docs/tutorials/index.md b/docs/tutorials/index.md
index bec4596..c4a5b3d 100644
--- a/docs/tutorials/index.md
+++ b/docs/tutorials/index.md
@@ -18,7 +18,6 @@ single_cell_eqtl_saigeqtl
 ehrdataset
 run_dataloader
 cell_level_ldsc_analysis
-ldsc_duncan_merge_anndatas
 magma_analysis_tutorial
 sclinker
 livi
diff --git a/docs/tutorials/ldsc_duncan_merge_anndatas.ipynb b/docs/tutorials/ldsc_duncan_merge_anndatas.ipynb
deleted file mode 100644
index cda2943..0000000
--- a/docs/tutorials/ldsc_duncan_merge_anndatas.ipynb
+++ /dev/null
@@ -1,286 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Merging AnnData objects for cell-level LDSC (Duncan lab dataset)\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "id": "3fdcf27c",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import anndata as ad\n",
-    "from pathlib import Path\n",
-    "from scipy.sparse import csr_array\n",
-    "import numpy as np\n",
-    "from tqdm.auto import tqdm"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 2,
-   "id": "00283b59",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import gc"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
-   "id": "22285917",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "obs_map = {\n",
-    "    \"Age\": \"age\",\n",
-    "    \"CellCycle\": \"cell_cycle\",\n",
-    "    \"CellID\": \"cell_id\",\n",
-    "    \"Chemistry\": \"chemistry\",\n",
-    "    \"Clusters\": \"clusters\",\n",
-    "    \"Donor\": \"donor\",\n",
-    "    \"NGenes\": \"n_genes\",\n",
-    "    \"ROIGroupCoarse\": \"roi_group_coarse\",\n",
-    "    \"ROIGroupFine\": \"roi_group_fine\",\n",
-    "    \"Roi\": \"roi\",\n",
-    "    \"SampleID\": \"sample_id\",\n",
-    "    \"Sex\": \"sex\",\n",
-    "    \"Subclusters\": \"subclusters\",\n",
-    "    \"Tissue\": \"tissue\",\n",
-    "    \"TotalUMI\": \"total_umi\",\n",
-    "}\n",
-    "\n",
-    "var_map = {\n",
-    "    \"Accession\": \"accession\",\n",
-    "    \"Gene\": \"gene\",\n",
-    "    \"Valid\": \"valid\",\n",
-    "}"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "id": "5b4b8160",
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "AnnData object with n_obs \u00d7 n_vars = 50000 \u00d7 59480\n",
-       "    obs: 'age', 'cell_cycle', 'chemistry', 'clusters', 'donor', 'n_genes', 'roi_group_coarse', 'roi_group_fine', 'roi', 'sample_id', 'sex', 'subclusters', 'tissue', 'total_umi', 'chunk'\n",
-       "    var: 'accession', 'gene', 'valid'"
-      ]
-     },
-     "execution_count": 4,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "adatas_paths = list(\n",
-    "    Path(\"/project/genomics/ayshan/ldsc_analysis/data_2/Single_cell_Siletti_Duncan/h5ad_chunks_new\").glob(\"*.h5ad\")\n",
-    ")\n",
-    "\n",
-    "\n",
-    "def read_h5ad(apath):\n",
-    "    adata = ad.read_h5ad(apath)\n",
-    "    adata.X = csr_array(adata.X.astype(np.uint16, copy=False))\n",
-    "    adata.obs = adata.obs.rename(columns=obs_map)\n",
-    "    adata.var = adata.var.rename(columns=var_map)\n",
-    "    adata.obs[\"chunk\"] = apath.stem\n",
-    "    adata.obs = adata.obs.set_index(\"cell_id\")\n",
-    "    return adata\n",
-    "\n",
-    "\n",
-    "adata = read_h5ad(adatas_paths[0])\n",
-    "adata.write_h5ad(\"/lustre/groups/shared/scgenetics/single_cell_siletti_duncan_single_chunk.h5ad\")\n",
-    "adata"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "id": "e0b2b8dd",
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "13f97890116b459394539a9ecec79ff1",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "  0%|          | 0/67 [00:00 **Where this connects to LDSC:** these same specificity scores were originally developed for stratified LD-score regression (S-LDSC), a separate method for partitioning trait heritability across annotations using GWAS summary statistics and a population reference panel. `cellink`'s `preprocess_for_sldsc` function computes the scores (hence the name), but the scores themselves — and everything in this notebook — are used independently of LDSC; MAGMA's gene-set and gene-property tests consume them directly. If you want the full S-LDSC heritability-partitioning workflow instead, see the companion notebook `cell_level_ldsc_analysis_updates.ipynb`.\n",
+    "> **Where this connects to LDSC:** these same specificity scores were originally developed for stratified LD-score regression (S-LDSC), a separate method for partitioning trait heritability across annotations using GWAS summary statistics and a population reference panel. `cellink`'s `preprocess_for_sldsc` function computes the scores (hence the name), but the scores themselves (and everything in this notebook) are used independently of LDSC; MAGMA's gene-set and gene-property tests consume them directly. If you want the full S-LDSC heritability-partitioning workflow instead, see the companion notebook `cell_level_ldsc_analysis_updates.ipynb`.\n",
     "\n",
-    "### GSA vs GPA — which one should I use?\n",
+    "### GSA vs GPA: which one should I use?\n",
     "\n",
     "Both tests start from the same specificity scores but ask the question differently:\n",
     "\n",
-    "- **GSA** thresholds each cell type's scores to a top-N% gene set (binary: in the set or not) and tests whether that set is enriched for GWAS signal. Simple to interpret, mirrors a standard binary enrichment test, and a good default — especially if you want results comparable to a published top-10%-style analysis.\n",
+    "- **GSA** thresholds each cell type's scores to a top-N% gene set (binary: in the set or not) and tests whether that set is enriched for GWAS signal. Simple to interpret, mirrors a standard binary enrichment test, and a good default, especially if you want results comparable to a published top-10%-style analysis.\n",
     "- **GPA** uses the full continuous score for every gene (no thresholding) and tests the linear relationship between score and GWAS gene-level association. It avoids the arbitrary top-N% cutoff and can be more powerful when specificity varies smoothly across genes rather than splitting cleanly into \"specific\" vs \"not.\"\n",
     "\n",
     "If you're unsure, start with GSA for its simplicity; switch to GPA if a hard top-N% cutoff feels like it's throwing away real signal in your data.\n",
@@ -160,10 +160,10 @@
     "# Path to the downloaded MAGMA binary (extracted by the cell above)\n",
     "MAGMA_BIN = os.path.join(\"magma_bin\", \"magma\")\n",
     "\n",
-    "# MAGMA gene location file — fetched from Ensembl and cached in ~/cellink_data/magma/\n",
+    "# MAGMA gene location file: fetched from Ensembl and cached in ~/cellink_data/magma/\n",
     "MAGMA_GENE_LOC = get_magma_gene_loc(genome_build=\"GRCh37\")\n",
     "\n",
-    "# 1000G PLINK files — downloaded and cached in ~/cellink_data/1000genomes_plink_EUR/\n",
+    "# 1000G PLINK files: downloaded and cached in ~/cellink_data/1000genomes_plink_EUR/\n",
     "_plink_path, _plink_prefix = get_1000genomes_plink_files(\n",
     "    config_path=str(Path(cellink.__file__).parent / \"resources\" / \"config\" / \"1000genomes.yaml\"),\n",
     "    population=\"EUR\",\n",
@@ -184,7 +184,7 @@
    "source": [
     "## Load and Prepare Data\n",
     "\n",
-    "We load the real OneK1K dataset and compute cell-type specificity scores as described above. (If you've already run the companion LDSC notebook, these preprocessing steps are identical — only the downstream analysis differs.)\n"
+    "We load the real OneK1K dataset and compute cell-type specificity scores as described above. (If you've already run the companion LDSC notebook, these preprocessing steps are identical; only the downstream analysis differs.)\n"
    ]
   },
   {
@@ -827,20 +827,20 @@
    "metadata": {},
    "source": [
     "---\n",
-    "## Part 1 — MAGMA \n",
+    "## Part 1: MAGMA \n",
     "\n",
     "MAGMA provides a complementary enrichment analysis that does not require LD score files. The pipeline has three steps:\n",
     "\n",
     "| Step | Function | Output |\n",
     "|---|---|---|\n",
-    "| I — Annotate | `run_magma_annotate` | `.genes.annot` |\n",
-    "| II — Gene analysis | `run_magma_gene_analysis` | `.genes.raw` |\n",
-    "| III-a — Gene-set analysis | `run_magma_gsa` | `.gsa.out` |\n",
-    "| III-b — Gene property analysis | `run_magma_gpa` | `.gsa.out` |\n",
+    "| I: Annotate | `run_magma_annotate` | `.genes.annot` |\n",
+    "| II: Gene analysis | `run_magma_gene_analysis` | `.genes.raw` |\n",
+    "| III-a: Gene-set analysis | `run_magma_gsa` | `.gsa.out` |\n",
+    "| III-b: Gene property analysis | `run_magma_gpa` | `.gsa.out` |\n",
     "\n",
     "Steps I and II only need to be run once per GWAS. Step III is run once per score set / method.\n",
     "\n",
-    "### Step I — Annotate SNPs to genes"
+    "### Step I: Annotate SNPs to genes"
    ]
   },
   {
@@ -992,7 +992,7 @@
    ],
    "source": [
     "# Derive MAGMA's SNP location file (SNP / CHR / BP) from real GWAS summary\n",
-    "# statistics — the same study (GCST004787, coronary artery disease) used in\n",
+    "# statistics, the same study (GCST004787, coronary artery disease) used in\n",
     "# the companion LDSC notebook's heritability/genetic-correlation sections.\n",
     "gwas_summary_statistic_path = get_gwas_catalog_study_summary_stats(\n",
     "    \"GCST004787\", genome_build=\"GRCh37\", return_path=True\n",
@@ -1000,7 +1000,7 @@
     "gwas_df = pd.read_csv(gwas_summary_statistic_path, sep=\"\\t\")\n",
     "gwas_df = gwas_df.dropna(subset=[\"variant_id\", \"chromosome\", \"base_pair_location\", \"p_value\"])\n",
     "# chromosome/base_pair_location load as float64 (NaN forces the column to\n",
-    "# float dtype) — cast to int before filtering/writing, otherwise \"22.0\" never\n",
+    "# float dtype), cast to int before filtering/writing, otherwise \"22.0\" never\n",
     "# matches \"22\" and MAGMA receives malformed CHR/BP columns.\n",
     "gwas_df[\"chromosome\"] = gwas_df[\"chromosome\"].astype(int)\n",
     "gwas_df[\"base_pair_location\"] = gwas_df[\"base_pair_location\"].astype(int)\n",
@@ -1061,7 +1061,7 @@
    "id": "cell-19",
    "metadata": {},
    "source": [
-    "### Step II — Gene-level association analysis\n",
+    "### Step II: Gene-level association analysis\n",
     "\n",
     "This step uses a PLINK LD reference panel to compute gene-level z-scores from the GWAS p-values. The output `.genes.raw` file is the input for both GSA and GPA."
    ]
@@ -1228,14 +1228,14 @@
    "metadata": {},
    "source": [
     "---\n",
-    "## Part 2 — MAGMA Gene-Set Analysis (GSA)\n",
+    "## Part 2: MAGMA Gene-Set Analysis (GSA)\n",
     "\n",
     "GSA tests whether the top-scoring genes in each cell type show stronger GWAS association than background. There are two paths to create the required GMT file:\n",
     "\n",
-    "* **From a scores DataFrame** — use `scores_to_gmt` (selects top-N% per cell type)\n",
-    "* **From existing LDSC `.GeneSet` files** — use `genesets_dir_to_entrez_gmt`\n",
+    "* **From a scores DataFrame**: use `scores_to_gmt` (selects top-N% per cell type)\n",
+    "* **From existing LDSC `.GeneSet` files**: use `genesets_dir_to_entrez_gmt`\n",
     "\n",
-    "### Path A — From continuous scores"
+    "### Path A: From continuous scores"
    ]
   },
   {
@@ -1318,7 +1318,7 @@
    "id": "cell-25",
    "metadata": {},
    "source": [
-    "### Path B — From existing binary LDSC `.GeneSet` files"
+    "### Path B: From existing binary LDSC `.GeneSet` files"
    ]
   },
   {
@@ -1692,9 +1692,9 @@
    "metadata": {},
    "source": [
     "---\n",
-    "## Part 3 — MAGMA Gene Property Analysis (GPA)\n",
+    "## Part 3: MAGMA Gene Property Analysis (GPA)\n",
     "\n",
-    "GPA tests the linear relationship between continuous per-gene scores and GWAS gene z-scores. Unlike GSA it uses all genes with scores — no top-N threshold. The covariate file contains one column per cell type.\n",
+    "GPA tests the linear relationship between continuous per-gene scores and GWAS gene z-scores. Unlike GSA it uses all genes with scores: no top-N threshold. The covariate file contains one column per cell type.\n",
     "\n",
     "### Build the covariate file"
    ]
@@ -2254,7 +2254,7 @@
    "source": [
     "### Univariate GPA\n",
     "\n",
-    "Each cell type is tested in an independent MAGMA call. Slower than joint mode but guarantees a result for every cell type — recommended when covariate columns are highly correlated (e.g. residual CV scores or negated scores)."
+    "Each cell type is tested in an independent MAGMA call. Slower than joint mode but guarantees a result for every cell type, recommended when covariate columns are highly correlated (e.g. residual CV scores or negated scores)."
    ]
   },
   {
@@ -4379,7 +4379,7 @@
    "metadata": {},
    "source": [
     "---\n",
-    "## Part 4 — LDSC → MAGMA: Full Workflows\n",
+    "## Part 4: LDSC → MAGMA: Full Workflows\n",
     "\n",
     "The two complete workflows from LDSC score outputs to MAGMA results.\n",
     "\n",
@@ -4647,7 +4647,7 @@
    "metadata": {},
    "source": [
     "---\n",
-    "## Part 5 — Tips and Common Patterns\n",
+    "## Part 5: Tips and Common Patterns\n",
     "\n",
     "### Previewing commands without running\n",
     "\n",
@@ -4865,9 +4865,9 @@
     "\n",
     "Key properties shared across all functions:\n",
     "\n",
-    "* **`run=False`** returns the command list without executing — useful for HPC job submission.\n",
-    "* **Gene ID handling** — `scores_to_gmt` and `scores_to_covar` accept a `gene_map` argument for symbol → ENSG translation; rows that cannot be mapped to ENSG are dropped.\n",
-    "* **Univariate GPA** — `run_magma_gpa(univariate=True)` tests each cell type independently, avoiding MAGMA's collinearity filter when scores are highly correlated."
+    "* **`run=False`** returns the command list without executing, useful for HPC job submission.\n",
+    "* **Gene ID handling**: `scores_to_gmt` and `scores_to_covar` accept a `gene_map` argument for symbol → ENSG translation; rows that cannot be mapped to ENSG are dropped.\n",
+    "* **Univariate GPA**: `run_magma_gpa(univariate=True)` tests each cell type independently, avoiding MAGMA's collinearity filter when scores are highly correlated."
    ]
   }
  ],
diff --git a/docs/tutorials/sclinker.ipynb b/docs/tutorials/sclinker.ipynb
index bd1f930..cd4c03e 100644
--- a/docs/tutorials/sclinker.ipynb
+++ b/docs/tutorials/sclinker.ipynb
@@ -213,7 +213,7 @@
    "source": [
     "### 1b: sc-linker enhancer-gene links and gene coordinates\n",
     "\n",
-    "These files are not part of the standard 1000G bundle — they come from the\n",
+    "These files are not part of the standard 1000G bundle; they come from the\n",
     "Broad sc-linker server and are specific to this pipeline."
    ]
   },
@@ -318,7 +318,7 @@
     "                print(f\"Using adata.var['{col}'] as gene_name\")\n",
     "                break\n",
     "        else:\n",
-    "            print(\"var_names are ENSG IDs — querying BioMart for HGNC mapping ...\")\n",
+    "            print(\"var_names are ENSG IDs: querying BioMart for HGNC mapping ...\")\n",
     "            from cellink.tl.external import get_gene_annotation\n",
     "\n",
     "            gene_coord = get_gene_annotation(gene_id_type=\"ensembl\")\n",
@@ -336,7 +336,7 @@
     "                _query_biomart_and_write_gene_coords(get_data_home())\n",
     "            print(\"  Done. HGNC names will be applied automatically in compute_celltype_programs.\")\n",
     "    else:\n",
-    "        print(\"adata.var['gene_name'] already present — will map ENSG→HGNC automatically\")\n",
+    "        print(\"adata.var['gene_name'] already present: will map ENSG→HGNC automatically\")\n",
     "\n",
     "print(f\"var_names sample: {adata.var_names[:3].tolist()}\")\n",
     "\n",
@@ -19270,7 +19270,7 @@
     "abc_road = results_df[results_df[\"strategy\"].str.startswith(\"ABC_Road\")].copy()\n",
     "\n",
     "if abc_road.empty:\n",
-    "    print(\"No ABC_Road results yet — run heritability first.\")\n",
+    "    print(\"No ABC_Road results yet: run heritability first.\")\n",
     "else:\n",
     "    programs = sorted(abc_road[\"program\"].unique())\n",
     "    traits = sorted(abc_road[\"trait\"].unique())\n",
diff --git a/docs/tutorials/scprs.ipynb b/docs/tutorials/scprs.ipynb
index 6d4ff11..2f44052 100644
--- a/docs/tutorials/scprs.ipynb
+++ b/docs/tutorials/scprs.ipynb
@@ -94,16 +94,16 @@
     "The code blocks below are **not executable cells**, they are reference code that you copy and run once (interactively in a Python session, or as a small CPU job on the cluster). After Phase 1 completes you return to this notebook for Phase 2.\n",
     "\n",
     "Running Phase 1 produces:\n",
-    "- `scprs_data/pbmc_atac_annotated.h5ad` — preprocessed ATAC object with peak coordinates\n",
-    "- `cluster_jobs/scprs_tutorial_array.sh` — SLURM array script (one PLINK task per cell)\n",
+    "- `scprs_data/pbmc_atac_annotated.h5ad`: preprocessed ATAC object with peak coordinates\n",
+    "- `cluster_jobs/scprs_tutorial_array.sh`: SLURM array script (one PLINK task per cell)\n",
     "\n",
     "Submit the array job with the printed `sbatch` command and wait for all tasks to finish.\n",
     "\n",
-    "> **Genome build — GWAS and reference panel must match.**\n",
+    "> **Genome build: GWAS and reference panel must match.**\n",
     ">\n",
     "> The critical constraint is that your GWAS summary statistics and your PLINK bfile must use the **same genome build**. Mixing builds causes position mismatches, silently near-zero PRS, and no error message.\n",
     ">\n",
-    "> You are not limited to 1000G or to the monocyte count GWAS used here — see *Adapting to your own analysis* at the end of this notebook."
+    "> You are not limited to 1000G or to the monocyte count GWAS used here; see *Adapting to your own analysis* at the end of this notebook."
    ]
   },
   {
@@ -138,7 +138,7 @@
    "source": [
     "### Phase 1b · RNA QC + cell type annotation\n",
     "\n",
-    "The QC thresholds and `cluster_to_celltype` mapping are specific to this PBMC dataset — adjust them based on your own dotplot. Cell type labels are **optional**: they are only used for the enrichment test in `get_disease_relevant_cells`. If you don't have matched RNA data, skip this step entirely and pass `cell_type_col=None` in Phase 2.\n",
+    "The QC thresholds and `cluster_to_celltype` mapping are specific to this PBMC dataset. Adjust them based on your own dotplot. Cell type labels are **optional**: they are only used for the enrichment test in `get_disease_relevant_cells`. If you don't have matched RNA data, skip this step entirely and pass `cell_type_col=None` in Phase 2.\n",
     "\n",
     "```python\n",
     "mdata = mu.read_10x_h5(str(H5_FILE))\n",
@@ -186,7 +186,7 @@
    "source": [
     "### Phase 1c · ATAC QC + LSI embedding\n",
     "\n",
-    "The `chrom`, `start`, `end` columns in `adata.var` are required by `write_slurm_array_job`. The parsing below assumes peak names in the format `chrN:start-end` (standard 10x output). Label transfer from RNA is optional — skip the last block if you have no cell type annotations.\n",
+    "The `chrom`, `start`, `end` columns in `adata.var` are required by `write_slurm_array_job`. The parsing below assumes peak names in the format `chrN:start-end` (standard 10x output). Label transfer from RNA is optional; skip the last block if you have no cell type annotations.\n",
     "\n",
     "```python\n",
     "atac = mdata[\"atac\"].copy()\n",
@@ -200,7 +200,7 @@
     "sc.pp.neighbors(atac, use_rep=\"X_lsi\", n_neighbors=25, n_pcs=20)\n",
     "sc.tl.umap(atac)\n",
     "\n",
-    "# Parse peak coordinates — required by write_slurm_array_job\n",
+    "# Parse peak coordinates: required by write_slurm_array_job\n",
     "parts = atac.var_names.str.replace(\":\", \"-\").str.split(\"-\")\n",
     "atac.var[\"chrom\"] = parts.str[0].str.replace(\"chr\", \"\", regex=False)\n",
     "atac.var[\"start\"] = parts.str[1].astype(int)\n",
@@ -241,13 +241,13 @@
     "                  names=[\"FID\", \"IID\", \"PAT\", \"MAT\", \"SEX\", \"PHENO\"])\n",
     "donor_ids = fam[\"IID\"].astype(str).tolist()\n",
     "np.random.seed(42)\n",
-    "labels = pd.Series(   # ← PLACEHOLDER — replace with real phenotypes\n",
+    "labels = pd.Series(   # ← PLACEHOLDER: replace with real phenotypes\n",
     "    np.random.choice([0, 1], size=len(donor_ids), p=[0.5, 0.5]),\n",
     "    index=donor_ids, name=\"label\",\n",
     ")\n",
     "\n",
     "# ── Tutorial-only: subsample for speed ───────────────────────────────────────\n",
-    "# Skip for real analyses — pass the full atac object to write_slurm_array_job.\n",
+    "# Skip for real analyses; pass the full atac object to write_slurm_array_job.\n",
     "atac_small = sc.read_h5ad(str(ATAC_H5AD))\n",
     "if atac_small.n_obs > N_CELLS_TUTORIAL:\n",
     "    sc.pp.subsample(atac_small, n_obs=N_CELLS_TUTORIAL, random_state=0)\n",
@@ -255,7 +255,7 @@
     "\n",
     "# ── SLURM array job ───────────────────────────────────────────────────────────\n",
     "# Adjust slurm_partition / slurm_qos for your cluster (check with `sinfo`).\n",
-    "# No SLURM? See Adapting section — the commands_file works with any scheduler.\n",
+    "# No SLURM? See Adapting section; the commands_file works with any scheduler.\n",
     "result = write_slurm_array_job(\n",
     "    adata=atac_small,           # real analysis: use atac (full object)\n",
     "    gwas_file=str(GWAS_FILE),\n",
@@ -313,7 +313,7 @@
     "if ATAC_SMALL_H5AD.exists():\n",
     "    atac_small = sc.read_h5ad(str(ATAC_SMALL_H5AD))\n",
     "else:\n",
-    "    print(f\"{ATAC_SMALL_H5AD.name} not found — reconstructing from {ATAC_H5AD.name}\")\n",
+    "    print(f\"{ATAC_SMALL_H5AD.name} not found, reconstructing from {ATAC_H5AD.name}\")\n",
     "    atac_small = sc.read_h5ad(str(ATAC_H5AD))\n",
     "    if atac_small.n_obs > N_CELLS_TUTORIAL:\n",
     "        sc.pp.subsample(atac_small, n_obs=N_CELLS_TUTORIAL, random_state=0)\n",
diff --git a/pyproject.toml b/pyproject.toml
index 085a4b8..b0e1226 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -35,7 +35,7 @@ dependencies = [
   "muon",
   "numexpr",
   "pandas",
-  "pandas-plink",
+  "pandas-plink>=2.3.0",
   "pandera",
   "polars",
   "regex",
diff --git a/src/cellink/_core/donordata.py b/src/cellink/_core/donordata.py
index a2bbfdd..75f6177 100644
--- a/src/cellink/_core/donordata.py
+++ b/src/cellink/_core/donordata.py
@@ -2,11 +2,13 @@
 
 import logging
 from collections.abc import Callable
-from dataclasses import dataclass
 
+import dask.array as da
 import h5py
+import numpy as np
 import pandas as pd
 import scanpy as sc
+import scipy.sparse as sp
 import zarr
 from anndata import AnnData
 from anndata.io import write_elem
@@ -25,7 +27,70 @@
 HIGHLIGHT_COLOR = "bold deep_pink2"
 
 
-@dataclass
+def _write_dense_array_zarr_chunked(f: zarr.Group, path: str, X, chunks: tuple[int, int] | None) -> None:
+    """Write one dense array (`X` or a dense layer) to `f[path]` with sane chunks."""
+    if isinstance(X, da.Array):
+        X = X.rechunk(chunks) if chunks is not None else X
+        da.to_zarr(X, path, overwrite=True)
+    else:
+        array_chunks = chunks if chunks is not None else tuple(min(dim, 4096) for dim in X.shape)
+        z = zarr.open_array(path, mode="w", shape=X.shape, chunks=array_chunks, dtype=X.dtype)
+        z[:] = np.asarray(X)
+    arr = zarr.open(path, mode="r+")
+    arr.attrs["encoding-type"] = "array"
+    arr.attrs["encoding-version"] = "0.2.0"
+
+
+def _write_anndata_zarr_dense_x_chunked(
+    f: zarr.Group, key: str, adata: AnnData, zarr_path: str, chunks: tuple[int, int] | None = None
+) -> None:
+    """Write one AnnData to a Zarr group with `X` and any dense `layers` chunked sanely.
+
+    `anndata.io.write_elem` picks its own chunk shape for a dense array
+    without regard to any chunking the input already has (confirmed: it
+    disregards a dask array's own `.chunks` entirely), which for a
+    genome-scale `X` (or a same-shaped dense layer, e.g. raw counts kept
+    alongside a normalized `X`) can pick a shape badly misaligned with how
+    the data is actually laid out on disk, making both the write and every
+    later read far slower than necessary. Writes everything except `X`
+    and any dense layers first (via a cheap, correctly-shaped zero-nnz
+    sparse placeholder for `X`, and by omitting dense layers entirely from
+    this first pass), then writes each dense array separately with an
+    explicit, sane `chunks`.
+
+    Parameters
+    ----------
+    chunks
+        Chunk shape applied to `X` and every dense layer. If `None`: a dask
+        array's own chunks (which `dask.array.to_zarr` preserves
+        automatically), or, for a plain array, capped at 4096 per axis.
+    """
+    X = adata.X
+    dense_layers = {k: v for k, v in adata.layers.items() if hasattr(v, "ndim") and not sp.issparse(v)}
+    sparse_layers = {k: v for k, v in adata.layers.items() if k not in dense_layers}
+
+    placeholder = sp.csr_matrix((adata.n_obs, adata.n_vars), dtype=X.dtype)
+    shell = AnnData(
+        X=placeholder,
+        obs=adata.obs,
+        var=adata.var,
+        obsm=dict(adata.obsm),
+        varm=dict(adata.varm),
+        obsp=dict(adata.obsp),
+        varp=dict(adata.varp),
+        layers=dict(sparse_layers),
+        uns=dict(adata.uns),
+    )
+    write_elem(f, key, shell)
+    del f[key]["X"]
+    _write_dense_array_zarr_chunked(f, f"{zarr_path}/{key}/X", X, chunks)
+
+    if dense_layers and "layers" not in f[key]:
+        f[key].create_group("layers")
+    for name, layer in dense_layers.items():
+        _write_dense_array_zarr_chunked(f, f"{zarr_path}/{key}/layers/{name}", layer, chunks)
+
+
 class DonorData:
     """Store and manage donor-related data with single-cell readouts.
 
@@ -56,6 +121,11 @@ def __init__(
             uns = {}
         if donor_id not in C.obs.columns:
             raise ValueError(f"'{donor_id}' not found in C.obs")
+        if isinstance(G, MuData) and G.obs.index.name is None:
+            for mod in G.mod.values():
+                if mod.obs.index.name == donor_id:
+                    G.obs.index.name = donor_id
+                    break
         if donor_id not in G.obs.columns and donor_id != G.obs.index.name:
             raise ValueError(f"'{donor_id}' must be in gdata.obs or set as index")
         if donor_id != G.obs.index.name:
@@ -103,17 +173,17 @@ def copy(self) -> DonorData:
             self._C = self._C.copy()
         return self
 
-    def _write_dd(self, f: h5py.File):
-        if isinstance(self.G, MuData):
-            g_group = f.create_group("G")
-            _write_h5mu(g_group, self.G)
-        else:
-            write_elem(f, "G", self.G)
-        if isinstance(self.C, MuData):
-            c_group = f.create_group("C")
-            _write_h5mu(c_group, self.C)
-        else:
-            write_elem(f, "C", self.C)
+    def _write_dd(self, f: h5py.File, zarr_path: str | None = None, x_chunks=None):
+        is_zarr = isinstance(f, zarr.Group)
+        for key, modality in (("G", self.G), ("C", self.C)):
+            if isinstance(modality, MuData):
+                group = f.create_group(key)
+                _write_h5mu(group, modality)
+            elif is_zarr and hasattr(modality.X, "ndim") and not sp.issparse(modality.X):
+                key_chunks = x_chunks.get(key) if isinstance(x_chunks, dict) else x_chunks
+                _write_anndata_zarr_dense_x_chunked(f, key, modality, zarr_path, chunks=key_chunks)
+            else:
+                write_elem(f, key, modality)
         f.attrs["encoding-type"] = "donordata"
 
         f.attrs["donor_id"] = self.donor_id
@@ -137,13 +207,22 @@ def write_h5_dd(self, path: str) -> None:
         with h5py.File(path, "w") as f:
             self._write_dd(f)
 
-    def write_zarr_dd(self, path: str) -> None:
+    def write_zarr_dd(
+        self, path: str, x_chunks: tuple[int, int] | dict[str, tuple[int, int]] | None = None
+    ) -> None:
         """Write the DonorData object to the specified file paths for both gene expression data (G) and cell-type data (C).
 
         Parameters
         ----------
         path : str | Path
             Path where the donor-data object should be saved.
+        x_chunks : tuple[int, int] | dict[str, tuple[int, int]], optional
+            Chunk shape for a dense `G.X`/`C.X`, e.g. `(4096, 4096)`. Pass a
+            dict (`{"G": ..., "C": ...}`) to set them independently. If not
+            given: a dask-backed `X` keeps its own existing chunks (via
+            `dask.array.to_zarr`); a plain array gets chunks capped at 4096
+            per axis. Only applies to a dense `X`; sparse/`MuData` are
+            unaffected.
 
         Example
         -------
@@ -153,7 +232,7 @@ def write_zarr_dd(self, path: str) -> None:
             if isinstance(m, MuData):
                 raise NotImplementedError("MuData not supported for zarr write")
         f = zarr.open(path, mode="w")
-        self._write_dd(f)
+        self._write_dd(f, zarr_path=str(path), x_chunks=x_chunks)
 
     def _ensure_extension(self, path: str, ext: str) -> str:
         """Ensure the given path ends with the desired extension."""
@@ -161,8 +240,8 @@ def _ensure_extension(self, path: str, ext: str) -> str:
             path += ext
         return path
 
-    def write_dd(self, path: str, dd: DonorData, fmt: str = None) -> None:
-        """Write the DonorData object to the specified file paths for both gene expression data (G) and cell-type data (C).
+    def write_dd(self, path: str, fmt: str = None) -> None:
+        """Write the DonorData object to the specified file path, format detected from the extension.
 
         Parameters
         ----------
@@ -183,10 +262,10 @@ def write_dd(self, path: str, dd: DonorData, fmt: str = None) -> None:
 
         if fmt == "h5":
             path = self._ensure_extension(path, ".dd.h5")
-            self.write_h5_dd(path, dd)
+            self.write_h5_dd(path)
         elif fmt == "zarr":
             path = self._ensure_extension(path, ".dd.zarr")
-            self.write_zarr_dd(path, dd)
+            self.write_zarr_dd(path)
         else:
             raise ValueError("Unknown format: use 'h5' or 'zarr'.")
 
@@ -430,8 +509,10 @@ def highlight_lines(lines, highlights):
 
     def __repr__(self) -> str:
         table = self.prep_repr()
-        Console().print(table)
-        return ""
+        console = Console()
+        with console.capture() as capture:
+            console.print(table)
+        return capture.get()
 
     def __str__(self) -> str:
         n_donors, n_donor_vars, n_cells, n_cell_vars = self.shape
diff --git a/src/cellink/cli/pgen.md b/src/cellink/cli/pgen.md
index aa387f4..66c3ead 100644
--- a/src/cellink/cli/pgen.md
+++ b/src/cellink/cli/pgen.md
@@ -23,7 +23,7 @@ cellink-pgen rare.pgen -o rare.zarr --sparse
 cellink-pgen rare.pgen common.pgen -o combined.zarr
 ```
 
-When multiple files are provided, variants are concatenated column-wise. All files must contain the **same samples in the same order** — obs is taken from the first file's `.psam`.
+When multiple files are provided, variants are concatenated column-wise. All files must contain the **same samples in the same order**: obs is taken from the first file's `.psam`.
 
 ## Loading the Output
 
@@ -41,7 +41,7 @@ X = adata.X.compute()
 
 The choice affects both write strategy and how `X` is stored on disk.
 
-**Dense** (default): Genotype blocks are streamed directly into a pre-allocated Zarr dataset. `X` is stored as a chunked `int8` array and returned as a Dask-backed array — reads remain lazy.
+**Dense** (default): Genotype blocks are streamed directly into a pre-allocated Zarr dataset. `X` is stored as a chunked `int8` array and returned as a Dask-backed array, so reads remain lazy.
 
 ```bash
 cellink-pgen common.pgen -o common.zarr
diff --git a/src/cellink/io/_pgen.py b/src/cellink/io/_pgen.py
index 8ecb8af..28e6f95 100644
--- a/src/cellink/io/_pgen.py
+++ b/src/cellink/io/_pgen.py
@@ -22,6 +22,40 @@
 logger = logging.getLogger(__name__)
 
 
+def _is_matrix_elem(elem_name: str) -> bool:
+    """True for ``X`` or any ``layers`` entry, the elements that can be a
+    dense, genome/transcriptome-scale array and so need lazy Dask loading.
+    """
+    return elem_name.endswith("/X") or elem_name.rsplit("/", 1)[0].endswith("/layers")
+
+
+def lazy_anndata_zarr_callback(func, elem_name: str, elem, iospec):
+    """``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
+    materializing it.
+    """
+    if iospec.encoding_type == "anndata" or elem_name.endswith("/"):
+        return ad.AnnData(
+            **{
+                k: read_dispatched(v, lazy_anndata_zarr_callback)
+                for k, v in dict(elem).items()
+                if not k.startswith("raw.")
+            }
+        )
+    elif _is_matrix_elem(elem_name) and iospec.encoding_type in (
+        "dataframe",
+        "csr_matrix",
+        "csc_matrix",
+        "awkward-array",
+    ):
+        return read_elem(elem)
+    elif _is_matrix_elem(elem_name) and iospec.encoding_type == "array":
+        return da.from_zarr(elem)
+    else:
+        return func(elem)
+
+
 def read_pgen_zarr(store: str | Path) -> ad.AnnData:
     """
     Lazily read an AnnData Zarr v3 store written by `stream_pgen_to_zarr`.
@@ -72,33 +106,15 @@ def read_pgen_zarr(store: str | Path) -> ad.AnnData:
     >>> X = adata.X.compute()
     """
     f = zarr.open(str(store), mode="r")
-
-    def callback(func, elem_name: str, elem, iospec):
-        if iospec.encoding_type == "anndata" or elem_name.endswith("/"):
-            return ad.AnnData(
-                **{k: read_dispatched(v, callback) for k, v in dict(elem).items() if not k.startswith("raw.")}
-            )
-        elif elem_name == "/X" and iospec.encoding_type in (
-            "dataframe",
-            "csr_matrix",
-            "csc_matrix",
-            "awkward-array",
-        ):
-            return read_elem(elem)
-        elif elem_name == "/X" and iospec.encoding_type == "array":
-            return da.from_zarr(elem)
-        else:
-            return func(elem)
-
-    return read_dispatched(f, callback=callback)
+    return read_dispatched(f, callback=lazy_anndata_zarr_callback)
 
 
 def _read_pvar(pvar_file: Path) -> pd.DataFrame:
-    """Read a PLINK2 .pvar file, correctly using its ``#CHROM POS ID REF ALT
-    ...`` header line for column names (mapped to cellink's canonical
-    ``chrom``/``pos``/``snp_id``/``a0``/``a1`` variant-annotation field
-    names, see :class:`cellink._core.data_fields.VAnn`) and skipping only
-    the preceding ``##`` metadata lines.
+    """Read a PLINK2 .pvar file, mapping its ``CHROM POS ID REF ALT ...``
+    columns to cellink's canonical ``chrom``/``pos``/``snp_id``/``a0``/``a1``
+    variant-annotation field names (see
+    :class:`cellink._core.data_fields.VAnn`), whether taken from a
+    ``#CHROM`` header line or, if absent, the standard PLINK2 column order.
     """
     n_meta = 0
     header_cols = None
@@ -112,18 +128,25 @@ def _read_pvar(pvar_file: Path) -> pd.DataFrame:
                 n_meta += 1
             break
 
+    standard_cols = ["CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"]
     if header_cols is None:
         pv = pd.read_csv(pvar_file, sep="\t", comment="#", header=None)
-        pv.columns = pv.columns.astype(str)
-        return pv
+        n = pv.shape[1]
+        header_cols = standard_cols[:n] if n <= len(standard_cols) else standard_cols + list(pv.columns[len(standard_cols) :])
+        pv.columns = header_cols
+    else:
+        pv = pd.read_csv(pvar_file, sep="\t", skiprows=n_meta, header=None, names=header_cols, dtype={"CHROM": str})
 
-    pv = pd.read_csv(pvar_file, sep="\t", skiprows=n_meta, header=None, names=header_cols, dtype={"CHROM": str})
     rename_map = {"CHROM": VAnn.chrom, "POS": VAnn.pos, "ID": VAnn.index, "REF": VAnn.a0, "ALT": VAnn.a1}
     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)
     if VAnn.chrom in pv.columns:
         pv[VAnn.chrom] = pv[VAnn.chrom].astype(str)
+    if VAnn.a0 in pv.columns:
+        pv[VAnn.a0] = pv[VAnn.a0].astype(str)
+    if VAnn.a1 in pv.columns:
+        pv[VAnn.a1] = pv[VAnn.a1].astype(str)
     return pv
 
 
diff --git a/src/cellink/io/_readwrite.py b/src/cellink/io/_readwrite.py
index 0144b7a..5bc3244 100644
--- a/src/cellink/io/_readwrite.py
+++ b/src/cellink/io/_readwrite.py
@@ -5,11 +5,13 @@
 from anndata._io.zarr import read_dataframe
 from anndata._types import StorageType
 from anndata.compat import _read_attr
+from anndata.experimental import read_dispatched
 from anndata.io import read_elem
 from mudata._core.io import _read_h5mu_mod
 from mudata._core.mudata import ModDict, MuData
 
 from cellink._core import DonorData
+from cellink.io._pgen import lazy_anndata_zarr_callback
 
 warnings.filterwarnings(
     "ignore",
@@ -72,9 +74,9 @@ def _read_mudata(group: StorageType, backed: bool = True) -> MuData:
     return mu
 
 
-def _read_dd(f: h5py.File) -> DonorData:
+def _read_dd(f: h5py.File, lazy: bool = False) -> DonorData:
     """
-    Internal function to read a DonorData object from an HDF5 file handle.
+    Internal function to read a DonorData object from an HDF5 or Zarr group.
 
     Reads donor-level genotype data (`G`) and cell-level expression data (`C`),
     reconstructs them as either AnnData or MuData objects, and collects additional
@@ -83,8 +85,13 @@ def _read_dd(f: h5py.File) -> DonorData:
 
     Parameters
     ----------
-    f : h5py.File
-        An open HDF5 file handle containing the DonorData object.
+    f : h5py.File or zarr.Group
+        An open HDF5 file handle or Zarr group containing the DonorData object.
+    lazy : bool, default=False
+        If True (Zarr only), read `G`/`C` via :func:`lazy_anndata_zarr_callback`
+        so a dense `X` stays Dask-backed instead of being materialized, which
+        lets callers select a subset (e.g. one cell type) before ever loading
+        the full array into memory.
 
     Returns
     -------
@@ -97,17 +104,21 @@ def _read_dd(f: h5py.File) -> DonorData:
     ValueError
         If the encoding type of `G` or `C` is not recognized.
     """
+
+    def _read_anndata(group):
+        return read_dispatched(group, callback=lazy_anndata_zarr_callback) if lazy else read_elem(group)
+
     if f["G"].attrs.get("encoding-type") == "MuData":
         G = _read_mudata(group=f["G"])
     elif f["G"].attrs.get("encoding-type") == "anndata":
-        G = read_elem(f["G"])
+        G = _read_anndata(f["G"])
     else:
         raise ValueError("Unknown encoding type for G")
 
     if f["C"].attrs.get("encoding-type") == "MuData":
         C = _read_mudata(group=f["C"])
     elif f["C"].attrs.get("encoding-type") == "anndata":
-        C = read_elem(f["C"])
+        C = _read_anndata(f["C"])
     else:
         raise ValueError("Unknown encoding type for C")
 
@@ -147,6 +158,11 @@ def read_zarr_dd(path: str) -> DonorData:
     """
     Read a DonorData object from a Zarr store on disk.
 
+    `G` and `C` are read lazily: a dense `X` stays Dask-backed rather than
+    being materialized on load, so a caller can select a subset (e.g. one
+    cell type via ``dd.sel(C_obs=...)``) and only that subset is ever read
+    off disk once it's actually computed/materialized.
+
     Parameters
     ----------
     path : str
@@ -158,7 +174,7 @@ def read_zarr_dd(path: str) -> DonorData:
         A DonorData object with genotype (`G`), cell expression (`C`), and metadata.
     """
     f = zarr.open(path, mode="r")
-    return _read_dd(f)
+    return _read_dd(f, lazy=True)
 
 
 def read_dd(path: str, fmt: str = None) -> DonorData:
diff --git a/src/cellink/io/_sgkit.py b/src/cellink/io/_sgkit.py
index a1bbe11..046867b 100644
--- a/src/cellink/io/_sgkit.py
+++ b/src/cellink/io/_sgkit.py
@@ -206,12 +206,12 @@ def from_sgkit_dataset(
     var_df.columns = var_df.columns.str.replace("variant_", "")
 
     if alleles_arr is not None:
-        var_df[VAnn.a0] = alleles_arr[:, 0]
+        var_df[VAnn.a0] = alleles_arr[:, 0].astype(str)
         if alleles_arr.shape[1] > 1:
-            var_df[VAnn.a1] = alleles_arr[:, 1]
+            var_df[VAnn.a1] = alleles_arr[:, 1].astype(str)
         if keep_multiallelic and alleles_arr.shape[1] > 2:
             for ai in range(2, alleles_arr.shape[1]):
-                var_df[f"a{ai}"] = alleles_arr[:, ai]
+                var_df[f"a{ai}"] = alleles_arr[:, ai].astype(str)
 
     if SgVars.contig_label in ds.data_vars:
         contigs = asarray(ds[SgVars.contig_label])
@@ -221,6 +221,7 @@ def from_sgkit_dataset(
                 var_df["chrom"] = pd.Index(vc).map(dict(enumerate(contigs)))
             except (KeyError, TypeError, ValueError):
                 var_df["chrom"] = vc
+            var_df["chrom"] = var_df["chrom"].astype(str)
         else:
             pass
 
@@ -263,7 +264,7 @@ def from_sgkit_dataset(
             for h in range(ploidy):
                 adata.layers[f"PHASE_{h}"] = phased_data[:, :, h].T
         else:
-            # 2D (variants, samples) — single phased flag per call
+            # 2D (variants, samples): single phased flag per call
             adata.layers["PHASE_0"] = phased_data.T
     else:
         adata.uns["has_phased_flag"] = False
diff --git a/src/cellink/pl/_ldsc_magma.py b/src/cellink/pl/_ldsc_magma.py
index 04a8b98..5d0d5f5 100644
--- a/src/cellink/pl/_ldsc_magma.py
+++ b/src/cellink/pl/_ldsc_magma.py
@@ -36,7 +36,7 @@ def enrichment_scatter(
     thresholds.  Significant cell types are listed in a legend outside the plot.
 
     Works identically for LDSC cell-type results (``Coefficient_P_value`` converted
-    to -log10) and MAGMA GSA / GPA results (``P`` column converted to -log10) —
+    to -log10) and MAGMA GSA / GPA results (``P`` column converted to -log10);
     the caller is responsible for building the DataFrame.
 
     Parameters
@@ -195,7 +195,7 @@ def method_comparison(
     """
     Compare -log10(p) values between two methods across shared cell types.
 
-    Produces a grid of scatter plots — one panel per shared cell type.  Within
+    Produces a grid of scatter plots, one panel per shared cell type.  Within
     each panel, each point represents a phenotype; its x-coordinate is the
     -log10(p) from ``df_a`` and its y-coordinate from ``df_b``.  A Spearman
     correlation coefficient and p-value are shown in each panel title.
diff --git a/src/cellink/resources/_gwas_prs_qtl.py b/src/cellink/resources/_gwas_prs_qtl.py
index 3bbeda4..4cd7e3c 100644
--- a/src/cellink/resources/_gwas_prs_qtl.py
+++ b/src/cellink/resources/_gwas_prs_qtl.py
@@ -27,16 +27,15 @@ def _normalize_build(genome_build: str) -> str:
 def _find_candidate_files(html: str) -> list[str]:
     """Find likely summary-stats filenames in an FTP directory listing.
 
-    Matches anything ending in .tsv, .txt, .zip, or .gz -- some
+    Matches anything ending in .tsv, .txt, .zip, or .gz, since some
     pre-harmonisation-era deposits ship a plain/zipped .txt instead of
-    .tsv.gz, or a bare ".gz" with no .tsv/.txt in the name at all (confirmed
-    directly: a sleep-duration GWAS whose only file is
-    "..._sumstats.txt.zip", and a major-depression GWAS whose only file is
-    literally "MDD2018_ex23andMe.gz"). Since .tsv.gz/.txt.gz already end in
-    ".gz", matching bare ".gz" covers all of those cases in one pattern.
-    Widening this far risks also matching README/LICENSE/CHANGELOG files
-    that live in the same directory (which the original .tsv.gz-only regex
-    never collided with) -- hence the explicit exclusion below.
+    .tsv.gz, or a bare ".gz" with no .tsv/.txt in the name at all (e.g. a
+    sleep-duration GWAS whose only file is "..._sumstats.txt.zip", or a
+    major-depression GWAS whose only file is literally
+    "MDD2018_ex23andMe.gz"). Since .tsv.gz/.txt.gz already end in ".gz",
+    matching bare ".gz" covers all of those cases in one pattern. Widening
+    this far also risks matching README/LICENSE/CHANGELOG files that live
+    in the same directory, hence the explicit exclusion below.
     """
     files = re.findall(r'href="([^"]*\.(?:tsv|txt|zip|gz))"', html)
     return [f for f in files if not _NON_DATA_FILENAME.search(f)]
@@ -171,7 +170,7 @@ def get_gwas_catalog_study_summary_stats(
         DataFrame containing the summary statistics, or Path to the downloaded file if return_path=True.
     """
     if translate_to_build and return_path:
-        raise ValueError("translate_to_build requires return_path=False — liftover operates on an in-memory DataFrame.")
+        raise ValueError("translate_to_build requires return_path=False, since liftover operates on an in-memory DataFrame.")
     study_meta = _fetch(f"{GWAS_API_BASE}/studies/{accession_id}", params=params, paginate=False)
 
     if "full_summary_stats" not in study_meta:
diff --git a/src/cellink/tl/__init__.py b/src/cellink/tl/__init__.py
index 681662b..3a07235 100644
--- a/src/cellink/tl/__init__.py
+++ b/src/cellink/tl/__init__.py
@@ -8,5 +8,6 @@
     run_snpeff,
     run_vep,
 )
+from ._coloc import DEFAULT_PRIOR_VAR_CC, coloc_abf, coloc_susie
 from ._rvat import beta_weighting, run_burden_test, run_skat_test
 from ._subset_region import subset_gene, subset_genomic_region
diff --git a/src/cellink/tl/_annotate_snps_ensembl_rest.py b/src/cellink/tl/_annotate_snps_ensembl_rest.py
index 658fd25..24dc19f 100644
--- a/src/cellink/tl/_annotate_snps_ensembl_rest.py
+++ b/src/cellink/tl/_annotate_snps_ensembl_rest.py
@@ -6,14 +6,11 @@
 
 
 def _split_snp_id(snp_str):
-    chrom = snp_str.split(":")[0]
-    pos = snp_str.split(":")[1].split("_")[0]
-    a0 = snp_str.split(":")[1].split("_")[1]
-    a1 = snp_str.split(":")[1].split("_")[2]
+    chrom, pos, a0, a1 = snp_str.split(":")
     return chrom, pos, a0, a1
 
 
-def get_snp_df(variant_codes, server="https://grch37.rest.ensembl.org/"):
+def get_snp_df(variant_codes, server="https://rest.ensembl.org/"):
     """
     Retrieve SNP (Single Nucleotide Polymorphism) information and overlap with genes from Ensembl.
 
@@ -25,9 +22,10 @@ def get_snp_df(variant_codes, server="https://grch37.rest.ensembl.org/"):
     Parameters
     ----------
     variant_codes : list of str
-        A list of SNP identifiers in the format of chromosome, position, and alleles (e.g., `1_55516888_T_C`).
+        A list of SNP identifiers as colon-separated chromosome, position, and alleles (e.g., `1:55516888:T:C`).
     server : str, optional
-        The URL of the Ensembl REST API server to query. Defaults to the GRCh37 Ensembl server.
+        The URL of the Ensembl REST API server to query. Defaults to the GRCh38 Ensembl server. Pass
+        `https://grch37.rest.ensembl.org/` explicitly for GRCh37-coordinate variant codes.
 
     Returns
     -------
@@ -45,12 +43,12 @@ def get_snp_df(variant_codes, server="https://grch37.rest.ensembl.org/"):
 
     Example
     -------
-    >>> variant_codes = ["1_55516888_T_C", "2_117900001_A_G"]
+    >>> variant_codes = ["1:55516888:T:C", "2:117900001:A:G"]
     >>> var_df, gene_df = get_snp_df(variant_codes)
     >>> var_df.head()
           snp_id  is_in_gene    genes  ...  clinical_significance
-    0  1_55516888_T_C        True    GENE1  ...          pathogenic
-    1  2_117900001_A_G       False     GENE2 ...                benign
+    0  1:55516888:T:C        True    GENE1  ...          pathogenic
+    1  2:117900001:A:G       False     GENE2 ...                benign
     >>> gene_df.head()
           biotype    start     end  ...       strand
     id
diff --git a/src/cellink/tl/_coloc.py b/src/cellink/tl/_coloc.py
new file mode 100644
index 0000000..682b534
--- /dev/null
+++ b/src/cellink/tl/_coloc.py
@@ -0,0 +1,240 @@
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+__all__ = ["coloc_abf", "coloc_susie", "DEFAULT_PRIOR_VAR", "DEFAULT_PRIOR_VAR_CC"]
+
+
+DEFAULT_P1 = 1e-4
+DEFAULT_P2 = 1e-4
+DEFAULT_P12 = 1e-5
+DEFAULT_PRIOR_VAR = 0.15**2
+DEFAULT_PRIOR_VAR_CC = 0.2**2
+
+
+def _log_abf(beta: np.ndarray, se: np.ndarray, prior_var: float) -> np.ndarray:
+    """Wakefield (2007) approximate log Bayes factor per SNP."""
+    z = beta / se
+    v = se**2
+    r = prior_var / (prior_var + v)
+    return 0.5 * (np.log(1 - r) + r * z**2)
+
+
+def _logsumexp(x: np.ndarray) -> float:
+    m = np.max(x)
+    return m + np.log(np.sum(np.exp(x - m)))
+
+
+def _log_h3(l1: np.ndarray, l2: np.ndarray, lsum1: float, lsum2: float) -> float:
+    log_all_pairs = lsum1 + lsum2
+    log_diag = _logsumexp(l1 + l2)
+    diff = np.exp(log_all_pairs - log_diag) - 1
+    if diff <= 0:
+        return -np.inf
+    return log_diag + np.log(diff)
+
+
+def coloc_abf(
+    beta1: np.ndarray | pd.Series,
+    se1: np.ndarray | pd.Series,
+    beta2: np.ndarray | pd.Series,
+    se2: np.ndarray | pd.Series,
+    p1: float = DEFAULT_P1,
+    p2: float = DEFAULT_P2,
+    p12: float = DEFAULT_P12,
+    prior_var: float = DEFAULT_PRIOR_VAR,
+    prior_var2: float | None = None,
+) -> dict[str, float]:
+    """Approximate Bayes factor colocalization between two association
+    studies sharing the same set of SNPs (Giambartolomei et al. 2014).
+
+    Tests, for one gene/region, whether a cis-eQTL signal and a GWAS
+    signal are more consistent with a single shared causal variant (H4)
+    than with two distinct causal variants (H3), association with only
+    one trait (H1/H2), or no association at all (H0).
+
+    Parameters
+    ----------
+    beta1
+        Effect sizes for study 1 (e.g. cis-eQTL), one per SNP.
+    se1
+        Standard errors for study 1, same order as ``beta1``.
+    beta2
+        Effect sizes for study 2 (e.g. GWAS), same SNPs and order as ``beta1``.
+    se2
+        Standard errors for study 2, same order as ``beta2``.
+    p1
+        Prior probability a given SNP is causal for study 1 only.
+    p2
+        Prior probability a given SNP is causal for study 2 only.
+    p12
+        Prior probability a given SNP is causal for both studies jointly
+        (i.e. the shared-causal-variant hypothesis, H4).
+    prior_var
+        Prior variance of study 1's effect size, on the scale ``beta1``/``se1``
+        are given in (the default assumes a standardized, sdY=1 quantitative
+        trait; rescale by ``sdY**2`` otherwise, matching R coloc's own
+        ``sdY`` convention).
+    prior_var2
+        Prior variance of study 2's effect size. Defaults to ``prior_var``
+        (both studies quantitative, sharing one prior, the common case for
+        two eQTL-like studies). Pass ``0.2**2`` when study 2 is a
+        case-control GWAS, matching R coloc's own default for ``type="cc"``.
+
+    Returns
+    -------
+    dict with keys ``PP0``-``PP4``, the posterior probability of each
+    hypothesis (summing to 1):
+
+    - PP0: no association with either trait
+    - PP1: association with study 1 only
+    - PP2: association with study 2 only
+    - PP3: both traits associated, but with distinct causal variants
+    - PP4: both traits associated, with the same causal variant
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from cellink.tl import coloc_abf
+    >>> rng = np.random.default_rng(0)
+    >>> n = 50
+    >>> beta1 = rng.normal(0, 0.02, n)
+    >>> se1 = np.full(n, 0.05)
+    >>> beta2 = rng.normal(0, 0.02, n)
+    >>> se2 = np.full(n, 0.05)
+    >>> beta1[0], se1[0] = 0.5, 0.04  # SNP 0: a real shared causal signal
+    >>> beta2[0], se2[0] = 0.4, 0.05
+    >>> pp = coloc_abf(beta1, se1, beta2, se2)
+    >>> bool(pp["PP4"] > 0.99)
+    True
+    """
+    beta1 = np.asarray(beta1, dtype=float)
+    se1 = np.asarray(se1, dtype=float)
+    beta2 = np.asarray(beta2, dtype=float)
+    se2 = np.asarray(se2, dtype=float)
+    if not (len(beta1) == len(se1) == len(beta2) == len(se2)):
+        raise ValueError("beta1, se1, beta2, se2 must all have the same length (one shared set of SNPs).")
+
+    l1 = _log_abf(beta1, se1, prior_var)
+    l2 = _log_abf(beta2, se2, prior_var2 if prior_var2 is not None else prior_var)
+    return _combine_log_abf(l1, l2, p1, p2, p12)
+
+
+def _combine_log_abf(l1: np.ndarray, l2: np.ndarray, p1: float, p2: float, p12: float) -> dict[str, float]:
+    """Standard coloc H0-H4 combination (Giambartolomei et al. 2014, matching
+    R coloc's own ``combine.abf``) given two arrays of per-SNP log Bayes
+    factors already on a shared set of SNPs, the common core of both
+    ``coloc_abf`` (l1/l2 from Wakefield ABF on beta/se) and ``coloc_susie``
+    (l1/l2 taken directly from SuSiE's own per-effect ``lbf_variable`` rows)."""
+    if len(l1) != len(l2):
+        raise ValueError("l1 and l2 must have the same length (one shared set of SNPs).")
+    lsum1 = _logsumexp(l1)
+    lsum2 = _logsumexp(l2)
+    lsum12 = _logsumexp(l1 + l2)
+
+    log_h = {
+        "PP0": 0.0,
+        "PP1": np.log(p1) + lsum1,
+        "PP2": np.log(p2) + lsum2,
+        "PP3": np.log(p1) + np.log(p2) + _log_h3(l1, l2, lsum1, lsum2),
+        "PP4": np.log(p12) + lsum12,
+    }
+    logs = np.array(list(log_h.values()))
+    m = np.max(logs)
+    post = np.exp(logs - m)
+    post = post / post.sum()
+    return dict(zip(log_h.keys(), post))
+
+
+def coloc_susie(
+    lbf1: np.ndarray | pd.DataFrame,
+    lbf2: np.ndarray | pd.DataFrame,
+    cs1_index: list[int] | np.ndarray,
+    cs2_index: list[int] | np.ndarray,
+    p1: float = DEFAULT_P1,
+    p2: float = DEFAULT_P2,
+    p12: float = DEFAULT_P12,
+    snp_ids: list[str] | None = None,
+) -> pd.DataFrame:
+    """Multiple-causal-variant colocalization against SuSiE credible sets
+    from both studies (Wallace 2021; R coloc's ``coloc.susie()`` /
+    ``coloc.bf_bf()``), rather than ``coloc_abf``'s single-causal-variant-
+    per-trait assumption.
+
+    Does not run SuSiE itself, takes each study's own already-fitted
+    SuSiE per-single-effect log Bayes factor matrix (``lbf_variable``, a
+    real, standard output field of ``susieR::susie()``/``susie_rss()``: one
+    row per single effect ``l``, one column per SNP) and, for every pair of
+    real credible sets (one from each study), runs the exact same H0-H4
+    combination ``coloc_abf`` uses, just on that effect's own per-SNP
+    log-BF row directly, rather than a Wakefield ABF recomputed from
+    beta/se. This is the identical algorithm R's ``coloc.bf_bf()``
+    implements, verified numerically against it directly.
+
+    Parameters
+    ----------
+    lbf1
+        Study 1's SuSiE fit's ``lbf_variable`` matrix, shape
+        ``(n_effects1, n_snps)``, same SNP order/columns as ``lbf2``.
+    lbf2
+        Study 2's SuSiE fit's ``lbf_variable`` matrix, shape
+        ``(n_effects2, n_snps)``.
+    cs1_index
+        Row indices into ``lbf1`` corresponding to study 1's real credible
+        sets (SuSiE's own ``fit$sets$cs_index``, i.e. the single effects
+        that survived SuSiE's own purity/coverage filtering, not every
+        row of ``lbf_variable``, which includes effects SuSiE fit but did
+        not resolve into a credible set).
+    cs2_index
+        Row indices into ``lbf2`` for study 2's real credible sets.
+    p1, p2, p12
+        Same priors as ``coloc_abf``.
+    snp_ids
+        Optional SNP labels (length ``n_snps``) for the ``hit1``/``hit2``
+        columns (each pair's most likely causal SNP per study). Defaults to
+        integer positions if not given.
+
+    Returns
+    -------
+    pd.DataFrame, one row per (study-1 signal, study-2 signal) pair, with
+    columns ``idx1``, ``idx2`` (the row indices from ``cs1_index``/
+    ``cs2_index``), ``nsnps``, ``hit1``, ``hit2``, and ``PP0``-``PP4``.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from cellink.tl import coloc_susie
+    >>> n = 30
+    >>> lbf1 = np.zeros((1, n)); lbf1[0, 5] = 20.0  # effect 0 points at SNP 5
+    >>> lbf2 = np.zeros((1, n)); lbf2[0, 5] = 15.0  # effect 0 also points at SNP 5
+    >>> res = coloc_susie(lbf1, lbf2, cs1_index=[0], cs2_index=[0])
+    >>> bool(res.loc[0, "PP4"] > 0.9)
+    True
+    """
+    lbf1 = np.asarray(lbf1, dtype=float)
+    lbf2 = np.asarray(lbf2, dtype=float)
+    cs1_index = np.asarray(cs1_index, dtype=int)
+    cs2_index = np.asarray(cs2_index, dtype=int)
+    if lbf1.shape[1] != lbf2.shape[1]:
+        raise ValueError("lbf1 and lbf2 must have the same number of columns (one shared set of SNPs).")
+    n_snps = lbf1.shape[1]
+    labels = list(snp_ids) if snp_ids is not None else list(range(n_snps))
+
+    rows = []
+    for idx1 in cs1_index:
+        l1 = lbf1[idx1]
+        for idx2 in cs2_index:
+            l2 = lbf2[idx2]
+            pp = _combine_log_abf(l1, l2, p1, p2, p12)
+            rows.append(
+                {
+                    "idx1": int(idx1),
+                    "idx2": int(idx2),
+                    "nsnps": n_snps,
+                    "hit1": labels[int(np.argmax(l1))],
+                    "hit2": labels[int(np.argmax(l2))],
+                    **pp,
+                }
+            )
+    return pd.DataFrame(rows)
diff --git a/src/cellink/tl/_runner.py b/src/cellink/tl/_runner.py
index dd94f77..8ea5392 100644
--- a/src/cellink/tl/_runner.py
+++ b/src/cellink/tl/_runner.py
@@ -176,7 +176,7 @@ def run_command(self, base_command: str, file_paths: list[str] = None, check: bo
           so the patched ``parse.py`` is always active.
         - **sandbox / rebuild**: ``singularity_image`` is updated in the runner
           config to point at the sandbox directory or rebuilt SIF, so
-          ``_build_container_command`` picks it up automatically — no special
+          ``_build_container_command`` picks it up automatically; no special
           handling needed here.
 
         Parameters
diff --git a/src/cellink/tl/external/__init__.py b/src/cellink/tl/external/__init__.py
index 3637e47..d82b1c6 100644
--- a/src/cellink/tl/external/__init__.py
+++ b/src/cellink/tl/external/__init__.py
@@ -79,7 +79,7 @@
 )
 from ._seismic import run_seismic
 from ._sldsc_utils import generate_gene_coord_file, generate_sldsc_genesets, get_magma_gene_loc, preprocess_for_sldsc
-from ._tensorqtl import read_tensorqtl_results, run_tensorqtl
+from ._tensorqtl import build_known_cis_eqtls_from_tensorqtl, read_tensorqtl_results, run_tensorqtl
 
 __all__ = [
     "read_jaxqtl_results",
@@ -88,6 +88,7 @@
     "calculate_pcs",
     "read_tensorqtl_results",
     "run_tensorqtl",
+    "build_known_cis_eqtls_from_tensorqtl",
     "run_mixmil",
     "run_scdrs",
     "run_seismic",
diff --git a/src/cellink/tl/external/_gsmap.py b/src/cellink/tl/external/_gsmap.py
index 2093a6f..677b64d 100644
--- a/src/cellink/tl/external/_gsmap.py
+++ b/src/cellink/tl/external/_gsmap.py
@@ -122,7 +122,7 @@ def format_gsmap_sumstats(
         # Writing the full DataFrame risks exposing multiple beta-like columns that gsMap
         # auto-detects, causing "Found 2 different BETA columns" errors. A common example:
         # GWAS Catalog harmonised files contain both 'hm_beta' and a plain 'beta' column
-        # with identical values — gsMap sees both as BETA and raises ValueError.
+        # with identical values, and gsMap sees both as BETA and raises ValueError.
         col_map = {
             snp: snp,
             a1: a1,
@@ -140,7 +140,7 @@ def format_gsmap_sumstats(
         cols_to_keep = [v for v in col_map.values() if v is not None and v in sumstats.columns]
         sumstats = sumstats[cols_to_keep].copy()
 
-        # Coerce beta-like columns to numeric — gsMap computes Z = sign(BETA) * sqrt(chi2.isf(P,1))
+        # Coerce beta-like columns to numeric, since gsMap computes Z = sign(BETA) * sqrt(chi2.isf(P,1))
         # and fails with a TypeError if BETA contains string values.
         for col in sumstats.columns:
             if col.lower() in {"beta", "hm_beta", "effect_size", "or", "odds_ratio"}:
@@ -261,7 +261,7 @@ def load_gsmap_results(
         "report_path": None,
     }
 
-    # Spatial LDSC — gsMap writes comma-separated .csv.gz
+    # Spatial LDSC: gsMap writes comma-separated .csv.gz
     ldsc_dir = sample_dir / "spatial_ldsc"
     if ldsc_dir.exists():
         ldsc_files = list(ldsc_dir.glob(f"*{trait_name}*.gz")) or list(ldsc_dir.glob("*.gz"))
diff --git a/src/cellink/tl/external/_ldsc2magma.py b/src/cellink/tl/external/_ldsc2magma.py
index 0181737..237acbf 100644
--- a/src/cellink/tl/external/_ldsc2magma.py
+++ b/src/cellink/tl/external/_ldsc2magma.py
@@ -309,9 +309,9 @@ def scores_to_gmt(
     gene_map
         Optional mapping from gene symbols to ENSG IDs.  Accepts:
 
-        * ``str`` or ``Path`` — path to a two-column TSV with headers
+        * ``str`` or ``Path``: path to a two-column TSV with headers
           ``gene_name`` and ``ensg_id``.
-        * ``pd.Series`` — index = gene symbol, values = ENSG ID.
+        * ``pd.Series``: index = gene symbol, values = ENSG ID.
 
         If provided, non-ENSG index entries are translated and rows that still
         do not look like ENSG IDs after mapping are dropped.
@@ -410,7 +410,7 @@ def scores_to_covar(
     out_file
         Output ``.covar`` file path.
     gene_map
-        Optional mapping from gene symbols to ENSG IDs — same formats as in
+        Optional mapping from gene symbols to ENSG IDs, using the same formats as in
         :func:`scores_to_gmt`.  Non-ENSG rows are dropped after mapping.
     negate
         If ``True``, multiply all scores by ``-1`` before writing.  Use this
@@ -477,7 +477,7 @@ def run_magma_annotate(
     **kwargs,
 ) -> "dict[str, Any]":
     """
-    Run MAGMA Step I — annotate SNPs to genes.
+    Run MAGMA Step I: annotate SNPs to genes.
 
     Maps each SNP to the gene(s) whose transcribed region (± ``window_kb`` kb)
     overlaps its position.  Creates a ``.genes.annot`` file consumed by
@@ -507,9 +507,9 @@ def run_magma_annotate(
     Returns
     -------
     dict
-        ``annot_file`` — path to ``.genes.annot``.
-        ``files_created`` — list of output paths (if ``run=True``).
-        ``command`` — command list (if ``run=False``).
+        ``annot_file``: path to ``.genes.annot``.
+        ``files_created``: list of output paths (if ``run=True``).
+        ``command``: command list (if ``run=False``).
 
     Examples
     --------
@@ -522,7 +522,7 @@ def run_magma_annotate(
 
     See Also
     --------
-    run_magma_gene_analysis : Step II — compute gene-level p-values.
+    run_magma_gene_analysis : Step II: compute gene-level p-values.
     """
     cmd = [magma_bin, "--annotate"]
     if window_kb:
@@ -552,7 +552,7 @@ def run_magma_gene_analysis(
     **kwargs,
 ) -> "dict[str, Any]":
     """
-    Run MAGMA Step II — gene-level association analysis.
+    Run MAGMA Step II: gene-level association analysis.
 
     Computes gene-level p-values and z-scores from GWAS SNP-level summary
     statistics, taking LD structure into account using a reference genotype
@@ -575,12 +575,12 @@ def run_magma_gene_analysis(
     n_samples
         Total GWAS sample size, used when every SNP shares the same N.
         Passed as ``N=`` in the ``--pval`` argument. Mutually
-        exclusive with ``ncol`` — MAGMA requires exactly one of the two.
+        exclusive with ``ncol``, since MAGMA requires exactly one of the two.
     ncol
         Name of a column in ``pval_file`` holding a per-SNP sample size
         (e.g. for a meta-analysis where N varies across SNPs). Passed as
         ``ncol=`` in the ``--pval`` argument. A column merely being
-        named ``"N"`` is **not** enough on its own — MAGMA always requires
+        named ``"N"`` is **not** enough on its own, because MAGMA always requires
         an explicit ``N=`` or ``ncol=`` modifier.
     magma_bin
         Path to the MAGMA binary.
@@ -592,9 +592,9 @@ def run_magma_gene_analysis(
     Returns
     -------
     dict
-        ``gene_results`` — path to ``.genes.raw``.
-        ``files_created`` — list of output paths (if ``run=True``).
-        ``command`` — command list (if ``run=False``).
+        ``gene_results``: path to ``.genes.raw``.
+        ``files_created``: list of output paths (if ``run=True``).
+        ``command``: command list (if ``run=False``).
 
     Examples
     --------
@@ -608,9 +608,9 @@ def run_magma_gene_analysis(
 
     See Also
     --------
-    run_magma_annotate : Step I — SNP-to-gene annotation.
-    run_magma_gsa : Step III — gene-set analysis.
-    run_magma_gpa : Step III — gene property analysis.
+    run_magma_annotate : Step I: SNP-to-gene annotation.
+    run_magma_gsa : Step III: gene-set analysis.
+    run_magma_gpa : Step III: gene property analysis.
     """
     cmd = [
         magma_bin,
@@ -655,7 +655,7 @@ def run_magma_gsa(
     **kwargs,
 ) -> "dict[str, Any]":
     """
-    Run MAGMA Step III — gene-set analysis (GSA).
+    Run MAGMA Step III: gene-set analysis (GSA).
 
     Tests whether genes in each set have higher GWAS association signals than
     background genes, using the gene-level results from
@@ -695,9 +695,9 @@ def run_magma_gsa(
     Returns
     -------
     dict
-        ``results_file`` — path to ``.gsa.out``.
-        ``files_created`` — list (if ``run=True``).
-        ``command`` — command list (if ``run=False``).
+        ``results_file``: path to ``.gsa.out``.
+        ``files_created``: list (if ``run=True``).
+        ``command``: command list (if ``run=False``).
 
     Examples
     --------
@@ -736,7 +736,7 @@ def run_magma_gpa(
     **kwargs,
 ) -> "dict[str, Any]":
     """
-    Run MAGMA Step III — gene property analysis (GPA).
+    Run MAGMA Step III: gene property analysis (GPA).
 
     Tests the linear association between continuous per-gene scores and GWAS
     gene-level z-scores.  Unlike GSA (which uses a top-N threshold), GPA uses
@@ -745,7 +745,7 @@ def run_magma_gpa(
     When ``univariate=False`` (default), all cell types are tested jointly in a
     single MAGMA call (``--gene-covar``).  This is efficient but MAGMA may drop
     highly collinear covariates.  Set ``univariate=True`` to test each cell
-    type independently — this is slower but always produces a result for every
+    type independently, which is slower but always produces a result for every
     cell type (matches ``run_magma_gpa_univariate.py`` in the pipeline).
 
     **LDSC → MAGMA GPA workflow**
@@ -777,9 +777,9 @@ def run_magma_gpa(
     Returns
     -------
     dict
-        ``results_file`` — path to ``.gsa.out``.
-        ``files_created`` — list (if ``run=True``).
-        ``command`` — command list (if ``run=False``, joint mode only).
+        ``results_file``: path to ``.gsa.out``.
+        ``files_created``: list (if ``run=True``).
+        ``command``: command list (if ``run=False``, joint mode only).
 
     Examples
     --------
@@ -820,7 +820,7 @@ def run_magma_gpa(
         results_file = f"{out_prefix}.gsa.out"
         return {"results_file": results_file, "files_created": [results_file, f"{out_prefix}.log"]}
 
-    # -- Univariate mode: one MAGMA call per cell-type covariate column -------
+    # Univariate mode: one MAGMA call per cell-type covariate column.
     covar = pd.read_csv(gene_covar, sep="\t", index_col=0)
     cell_types = covar.columns.tolist()
     logger.info("GPA univariate: %d cell types", len(cell_types))
diff --git a/src/cellink/tl/external/_livi.py b/src/cellink/tl/external/_livi.py
index 25365f2..f0241e4 100644
--- a/src/cellink/tl/external/_livi.py
+++ b/src/cellink/tl/external/_livi.py
@@ -23,7 +23,7 @@ class LIVIRunner:
         ``configs/`` directories).
     execution_mode : {"python_api", "subprocess"}
         ``"python_api"`` (default) imports LIVI directly into the current
-        process — no Hydra overhead, returns Python objects.
+        process, avoiding Hydra overhead and returning Python objects directly.
         ``"subprocess"`` runs LIVI's ``src/train.py`` as a child process via
         the Hydra CLI, which is better suited for isolated HPC job submission.
     python_executable : str
@@ -651,7 +651,7 @@ def train_livi(
         ``pytorch_lightning.Trainer``.
     enable_checkpointing : bool
         If *False*, skip the ``ModelCheckpoint`` callback entirely (and the
-        return value becomes *None*) -- useful for short benchmark/smoke runs
+        return value becomes *None*), useful for short benchmark/smoke runs
         where you don't want checkpoint I/O.
     enable_logger : bool
         Passed through to ``pytorch_lightning.Trainer`` as ``logger=``; set
@@ -859,19 +859,19 @@ def infer_livi(
         corresponding model component was not trained):
 
         ``"cell_state_latent"``
-            (n_cells × z_dim) — cell-state latent factors per cell.
+            (n_cells × z_dim): cell-state latent factors per cell.
         ``"cell_state_decoder"``
-            (n_genes × z_dim) — gene loadings of the cell-state decoder.
+            (n_genes × z_dim): gene loadings of the cell-state decoder.
         ``"D_embedding"``
-            (n_donors × n_DxC) — donor × cell-state interaction embeddings.
+            (n_donors × n_DxC): donor × cell-state interaction embeddings.
         ``"DxC_decoder"``
-            (n_genes × n_DxC) — gene loadings of the DxC decoder.
+            (n_genes × n_DxC): gene loadings of the DxC decoder.
         ``"V_embedding"``
-            (n_donors × n_persistent) — persistent donor factor embeddings.
+            (n_donors × n_persistent): persistent donor factor embeddings.
         ``"V_decoder"``
-            (n_genes × n_persistent) — gene loadings of the V decoder.
+            (n_genes × n_persistent): gene loadings of the V decoder.
         ``"assignment_matrix"``
-            (z_dim × n_DxC) — assignment matrix *A* mapping cell-state
+            (z_dim × n_DxC): assignment matrix *A* mapping cell-state
             factors to DxC factors.
 
     Examples
@@ -901,7 +901,7 @@ def infer_livi(
     )
     model.eval()
 
-    # Factorize donors — order must match training factorisation
+    # Factorize donors (order must match training factorisation)
     _, y_index = pd.factorize(adata.obs[individual_col], sort=False, use_na_sentinel=False)
 
     dataset = LIVIDataset(
@@ -944,11 +944,11 @@ def infer_livi(
 
     results: dict[str, pd.DataFrame] = {"cell_state_latent": cell_state_latent}
 
-    # Cell-state decoder weights — shape (n_genes, z_dim)
+    # Cell-state decoder weights, shape (n_genes, z_dim)
     cs_dec_np = model.decoder.mean[0].weight.detach().cpu().numpy()
     results["cell_state_decoder"] = pd.DataFrame(cs_dec_np, index=adata.var.index, columns=z_factor_cols)
 
-    # Donor embeddings — access weight tensors directly (avoids n_cells overhead)
+    # Donor embeddings: access weight tensors directly (avoids n_cells overhead)
     n_unique_donors = len(y_index)
 
     if model.n_DxC_factors != 0:
@@ -1063,7 +1063,7 @@ def run_livi_association_testing(
 
     Examples
     --------
-    Pass DonorData — genotype matrix, kinship, and gPCs are auto-extracted:
+    Passing DonorData auto-extracts the genotype matrix, kinship, and gPCs:
 
     >>> assoc = cl.tl.external.run_livi_association_testing(
     ...     results,
@@ -1129,7 +1129,7 @@ def run_livi_association_testing(
         genotype_pcs = pd.read_csv(genotype_pcs, index_col=0, sep=sep)
         genotype_pcs.index = genotype_pcs.index.astype(str)
 
-    # livi_testing always does covariates["intercept"] = 1.0 — it must be a DataFrame.
+    # livi_testing always does covariates["intercept"] = 1.0, so it must be a DataFrame.
     # Initialise with the donor index from D_context (or V_persistent) if not supplied.
     if covariates is None:
         ref = D_context if D_context is not None else V_persistent
diff --git a/src/cellink/tl/external/_livi_annbatch.py b/src/cellink/tl/external/_livi_annbatch.py
index 9e78baa..2f0edbc 100644
--- a/src/cellink/tl/external/_livi_annbatch.py
+++ b/src/cellink/tl/external/_livi_annbatch.py
@@ -29,7 +29,7 @@ def build_annbatch_collection(
     Parameters
     ----------
     c_group
-        The cell-level AnnData store to stream from -- a `zarr.Group` /
+        The cell-level AnnData store to stream from: a `zarr.Group` /
         `h5py.Group` (e.g. the "C" group of a DonorData `.dd.zarr`/`.dd.h5`
         cache: ``zarr.open(path, "r")["C"]``) or a path to a standalone
         AnnData file.
@@ -178,7 +178,7 @@ def _presence(gdata, snp_idx, donor_categories):
         if n_missing:
             logger.warning(
                 "%d missing genotype calls (%.2f%%) among the selected cis SNPs; treating them as "
-                "allele-absent (0) for this binary presence feature -- NaN > 0 is always False in "
+                "allele-absent (0) for this binary presence feature, since NaN > 0 is always False in "
                 "numpy, so this was previously happening silently.",
                 n_missing,
                 100 * n_missing / raw.size,
@@ -406,20 +406,20 @@ def train_livi_annbatch(
     `DatasetCollection` (built from `c_group` if `collection_path` doesn't
     exist yet) instead of being held fully in memory. Genotype (donors x
     SNPs) is still loaded into memory and looked up per cell, exactly like
-    `train_livi`'s `eqtl_genotypes` -- it's orders of magnitude smaller than
+    `train_livi`'s `eqtl_genotypes`, since it's orders of magnitude smaller than
     cell expression.
 
     Parameters
     ----------
     c_group
-        Cell-level AnnData store to stream from -- e.g. the "C" group of a
+        Cell-level AnnData store to stream from, e.g. the "C" group of a
         DonorData `.dd.zarr`/`.dd.h5` cache (``zarr.open(path, "r")["C"]``).
         Ignored if `collection_path` already has a built collection.
     gdata : AnnData
         Genotype (donors x SNPs) AnnData, e.g. from :func:`read_g_from_dd_store`.
     output_dir
         Unused for checkpointing currently (annbatch path doesn't checkpoint
-        by default -- see Notes); kept for signature parity with `train_livi`.
+        by default; see Notes), kept for signature parity with `train_livi`.
     collection_path
         Path to the (possibly to-be-built) annbatch `DatasetCollection`.
     donor_key
@@ -429,7 +429,7 @@ def train_livi_annbatch(
     known_cis_eqtls
         Real SNPs x genes 0/1 mapping (e.g. from genomic-distance windows or a
         curated association list). If `None`, a synthetic random mapping is
-        generated from `cis_snps`/`target_genes`/`cis_genes_per_snp` -- fine
+        generated from `cis_snps`/`target_genes`/`cis_genes_per_snp`; fine
         for smoke testing, not biologically meaningful.
     cis_snps, target_genes, cis_genes_per_snp
         See `CisGenotype`; only used when `known_cis_eqtls` is `None`.
@@ -466,7 +466,7 @@ def train_livi_annbatch(
 
     Notes
     -----
-    No `ModelCheckpoint` is added by default (unlike `train_livi`) -- pass one
+    No `ModelCheckpoint` is added by default (unlike `train_livi`); pass one
     via `callbacks` if you need checkpointing.
     """
     if runner is None:
diff --git a/src/cellink/tl/external/_magma.py b/src/cellink/tl/external/_magma.py
index 9737ae7..8424a95 100644
--- a/src/cellink/tl/external/_magma.py
+++ b/src/cellink/tl/external/_magma.py
@@ -313,7 +313,7 @@ def prepare_magma_inputs(
         - 'dd_genotypes'   : export genotypes from a DonorData object (requires `dd`)
         - 'reference_panel': download a 1000G panel (requires `reference_panel`)
         - 'external'       : use pre-built PLINK files you already have (requires `external_ld_prefix`)
-        - None             : skip LD reference — you must supply it yourself to run_magma_gene_analysis
+        - None             : skip LD reference, you must supply it yourself to run_magma_gene_analysis
     dd : DonorData, optional
         Required only when ld_source='dd_genotypes'.
     reference_panel : str, optional
@@ -597,11 +597,11 @@ def run_magma_pipeline(
     n_samples : int
         GWAS sample size.
     ld_source : str or None
-        LD reference strategy — see prepare_magma_inputs for full docs.
+        LD reference strategy; see prepare_magma_inputs for full docs.
         - 'dd_genotypes'   : use genotypes from DonorData (pass `dd`)
         - 'reference_panel': download a 1000G panel (pass `reference_panel`)
         - 'external'       : use existing PLINK files (pass `external_ld_prefix`)
-        - None             : raises an error — LD reference is required for gene analysis
+        - None             : raises an error, LD reference is required for gene analysis
     dd : DonorData, optional
         Required when ld_source='dd_genotypes'.
     reference_panel : str, optional
diff --git a/src/cellink/tl/external/_sclinker.py b/src/cellink/tl/external/_sclinker.py
index f12d37b..83fd687 100644
--- a/src/cellink/tl/external/_sclinker.py
+++ b/src/cellink/tl/external/_sclinker.py
@@ -288,9 +288,9 @@ def compute_nmf_programs(
     device : str, default ``"cuda"``
         Device for torchnmf backend: ``"cuda"`` or ``"cpu"``.
 
-        - ``"cuda"`` — uses GPU if available, raises a clear warning if CUDA is
+        - ``"cuda"``: uses GPU if available, raises a clear warning if CUDA is
           not found and falls back to CPU.
-        - ``"cpu"``  — forces CPU even if a GPU is present.
+        - ``"cpu"``: forces CPU even if a GPU is present.
 
         If ``torchnmf`` is not installed at all, cellink logs an install hint
         and falls back to sklearn NMF (which is slower but always available).
@@ -342,7 +342,7 @@ def compute_nmf_programs(
         from torchnmf.nmf import NMF as TorchNMF
     except ImportError:
         logger.warning(
-            "torchnmf is not installed — falling back to sklearn NMF, which is "
+            "torchnmf is not installed, falling back to sklearn NMF, which is "
             "significantly slower on large matrices (>50k cells).\n"
             "Install the faster backend with:\n"
             "  pip install torchnmf"
@@ -380,8 +380,8 @@ def compute_nmf_programs(
         )
         model = NMF(
             n_components=n_components,
-            init="nndsvda",  # truncated-SVD warm start — ~5-10x faster than "random"
-            solver="cd",  # coordinate descent — faster than multiplicative update
+            init="nndsvda",  # truncated-SVD warm start, ~5-10x faster than "random"
+            solver="cd",  # coordinate descent, faster than multiplicative update
             max_iter=500,
             tol=1e-4,
             random_state=random_state,
diff --git a/src/cellink/tl/external/_sclinker_utils.py b/src/cellink/tl/external/_sclinker_utils.py
index e339ab7..24dcee4 100644
--- a/src/cellink/tl/external/_sclinker_utils.py
+++ b/src/cellink/tl/external/_sclinker_utils.py
@@ -72,11 +72,11 @@ def download_sclinker_enhancer_links(
     out_dir
         Directory to download files into.
     tissue
-        Accepted for API compatibility but ignored — the combined files cover
+        Accepted for API compatibility but ignored; the combined files cover
         all tissues. Pass ``tissue`` to :func:`load_roadmap_links` /
         :func:`load_abc_links` to filter after downloading.
     chromosomes
-        Unused — kept for API consistency.
+        Unused, kept for API consistency.
 
     Returns
     -------
@@ -266,7 +266,7 @@ def download_sclinker_references(
     Download sc-linker reference files.
 
     **For 1000G files (PLINK, LD scores, weights, frq, HapMap3), use the
-    cellink resource helpers instead** — they download from Zenodo
+    cellink resource helpers instead**, since they download from Zenodo
     (https://zenodo.org/records/10515792) which is reliable and fast:
 
     .. code-block:: python
@@ -293,7 +293,7 @@ def download_sclinker_references(
     download_roadmap, download_abc, download_gene_coords
         Which enhancer-gene link files to download (default: all True).
     download_bims, download_frq, download_weights, download_hapmap3, download_baseline
-        Kept for API compatibility but ignored — use the cellink resource
+        Kept for API compatibility but ignored; use the cellink resource
         helpers listed above instead.
     chromosomes
         Ignored (kept for API compatibility).
@@ -595,7 +595,7 @@ def get_gene_annotation(
 
         ``"ensembl"`` *(default)*
             Ensembl stable IDs (e.g. ``ENSG00000099338``).  Use this when
-            your AnnData ``var_names`` are ENSG IDs — the typical case for
+            your AnnData ``var_names`` are ENSG IDs, the typical case for
             sc-linker gene programs derived from standard scRNA-seq pipelines.
 
         ``"hgnc"``
@@ -656,7 +656,7 @@ def load_gene_annotation(
         Explicit path to an existing coord file. Passed to
         :func:`get_gene_annotation` as ``path``.
     gene_id_type : ``"ensembl"`` | ``"hgnc"``
-        Which identifier is in the ``GENE`` column — must match your data.
+        Which identifier is in the ``GENE`` column; must match your data.
         Default is ``"ensembl"`` because sc-linker AnnData objects typically
         have ENSG ``var_names``.
     data_home
@@ -1353,9 +1353,9 @@ def load_sclinker_heritability_results(
     if not rows:
         logger.warning(
             "No results parsed. Possible causes:\n"
-            "  1. No .results files alongside the .log files — check that LDSC ran with "
+            "  1. No .results files alongside the .log files: check that LDSC ran with "
             "--overlap-annot and that annotation_prefixes= was passed to run_sclinker_heritability().\n"
-            "  2. LDSC exited before writing .results — open a .log file and check for errors."
+            "  2. LDSC exited before writing .results: open a .log file and check for errors."
         )
         return pd.DataFrame()
 
@@ -1446,11 +1446,11 @@ def _get_parse_py_path(runner) -> str | None:
 
     Resolution order:
 
-    1. ``runner.parse_script`` / ``runner.config["parse_script"]`` — the
+    1. ``runner.parse_script`` / ``runner.config["parse_script"]``: the
        canonical field on ``LDSCRunner``.  Always set this explicitly for
        Docker and Singularity; the path inside the container is typically
        ``/ldsc/ldscore/parse.py``.
-    2. Auto-discovery via PATH — find ``ldsc.py`` on the host PATH and infer
+    2. Auto-discovery via PATH: find ``ldsc.py`` on the host PATH and infer
        the sibling ``ldscore/parse.py``.  Works for local installs; will not
        work inside containers.
 
@@ -1529,15 +1529,15 @@ def _write_parse_py_via_runner(runner, patched_source: str) -> bool:
         Overwrites ``parse.py`` on disk; backs up original as ``parse.py.bak``.
 
     **docker**
-        ``docker cp`` + ``docker commit`` — the patched file is baked
+        ``docker cp`` + ``docker commit``: the patched file is baked
         permanently into the image. Idempotent and survives restarts.
 
     **singularity**
-        SIF images are read-only squashfs archives — they cannot be edited
+        SIF images are read-only squashfs archives, so they cannot be edited
         in-place.  cellink supports three strategies, chosen via
         ``singularity_patch_strategy`` in the runner config:
 
-        ``"overlay"`` *(default — HPC-friendly, no root needed)*
+        ``"overlay"`` *(default, HPC-friendly, no root needed)*
             Creates a persistent ext3 overlay image alongside the SIF.
             On every ``ldsc.py`` call cellink appends
             ``--overlay `` so the patch is always active.
@@ -1793,7 +1793,7 @@ def check_and_patch_ldsc_parse_bug(runner) -> dict:
 
     Works for **local**, **Docker** and **Singularity** execution modes.
 
-    **Configuration** — add ``parse_script`` to your runner config so cellink
+    **Configuration**: add ``parse_script`` to your runner config so cellink
     knows exactly where ``parse.py`` lives, especially for containerised setups:
 
     .. code-block:: python
@@ -1812,12 +1812,12 @@ def check_and_patch_ldsc_parse_bug(runner) -> dict:
 
     Behaviour by mode
     -----------------
-    - **local** — patches ``parse.py`` on disk and backs up the original as
+    - **local**: patches ``parse.py`` on disk and backs up the original as
       ``parse.py.bak``.
-    - **docker** — patches the file inside the image via ``docker cp`` +
+    - **docker**: patches the file inside the image via ``docker cp`` +
       ``docker commit``.  The image is updated in-place and the patch survives
       container restarts.
-    - **singularity** — SIF images are read-only, so cellink writes the patched
+    - **singularity**: SIF images are read-only, so cellink writes the patched
       ``parse.py`` to a host-side shadow directory (default:
       ``~/.cellink/ldsc_patch/ldscore/``) and injects
       ``PYTHONPATH=:$PYTHONPATH`` into every subsequent ``ldsc.py``
@@ -1832,11 +1832,11 @@ def check_and_patch_ldsc_parse_bug(runner) -> dict:
     Returns
     -------
     dict with keys:
-        ``"status"``     — ``"already_patched"``, ``"patched"``,
+        ``"status"``     : ``"already_patched"``, ``"patched"``,
                            ``"patch_failed"`` or ``"not_found"``
-        ``"mode"``       — ``"local"``, ``"docker"`` or ``"singularity"``
-        ``"parse_path"`` — resolved path to ``parse.py`` (or None)
-        ``"detail"``     — human-readable explanation
+        ``"mode"``       : ``"local"``, ``"docker"`` or ``"singularity"``
+        ``"parse_path"`` : resolved path to ``parse.py`` (or None)
+        ``"detail"``     : human-readable explanation
     """
     mode = runner.config.get("execution_mode", "local")
     parse_path = _get_parse_py_path(runner)
@@ -1859,7 +1859,7 @@ def check_and_patch_ldsc_parse_bug(runner) -> dict:
             "status": "already_patched",
             "mode": mode,
             "parse_path": parse_path,
-            "detail": "parse.py already patched (chr_ld[0].columns reindex present) — no action needed.",
+            "detail": "parse.py already patched (chr_ld[0].columns reindex present), no action needed.",
         }
 
     patched = source
@@ -1953,7 +1953,7 @@ def _merge_bedgraph_bedtools(bg: pd.DataFrame) -> pd.DataFrame:
 
 
 def _merge_bedgraph_python(bg: pd.DataFrame) -> pd.DataFrame:
-    """Interval merge with score summing — numpy-based, no iterrows."""
+    """Interval merge with score summing, numpy-based, no iterrows."""
     out_chrs, out_starts, out_ends, out_scores = [], [], [], []
     for chrom, grp in bg.groupby("chr", sort=False):
         starts = grp["start"].to_numpy(dtype=np.int64)
@@ -2148,7 +2148,7 @@ def _parse_ldsc_log(log_file: Path) -> dict | None:
     and runtime metadata to ``.log``.  This function reads both.
 
     Always extracted from ``.log`` (if present):
-        ``h2_obs``, ``h2_obs_se`` — observed-scale heritability.
+        ``h2_obs``, ``h2_obs_se``: observed-scale heritability.
 
     Extracted from ``.results`` (when ``--overlap-annot`` succeeded):
         ``Enrichment``, ``Enrichment_std_error``, ``Enrichment_z_score``,
diff --git a/src/cellink/tl/external/_scooby.py b/src/cellink/tl/external/_scooby.py
index f8b5fb5..863454a 100644
--- a/src/cellink/tl/external/_scooby.py
+++ b/src/cellink/tl/external/_scooby.py
@@ -1,6 +1,8 @@
 from __future__ import annotations
 
+import json
 import logging
+import re
 from pathlib import Path
 from typing import Any, Literal
 
@@ -234,12 +236,35 @@ def build_scooby_embedding_scpoli(
         _ad._cellink_read_shim_applied = True
     from scarches.models.scpoli import scPoli
 
+    import gc
+
+    import scipy.sparse
+
     a = adata
     if use_hvg:
+        assert scipy.sparse.issparse(a.X), (
+            f"expected a.X to stay sparse into the HVG step, got {type(a.X)}, "
+            "a dense copy here at full-cohort scale (7M+ cells) would itself "
+            "exhaust available memory long before any real peak."
+        )
+        sc.pp.filter_genes(a, min_cells=10)
+        gc.collect()
         sc.pp.highly_variable_genes(a, flavor="seurat_v3", n_top_genes=n_top_genes, batch_key=condition_key)
+        gc.collect()
+        a.X = a.X.tocsc()
         a = a[:, a.var["highly_variable"]].copy()
+        gc.collect()
+        adata.X = scipy.sparse.csr_matrix((adata.n_obs, adata.n_vars), dtype=adata.X.dtype)
+        gc.collect()
 
-    if checkpoint_dir is not None and Path(checkpoint_dir).exists():
+    if scipy.sparse.issparse(a.X):
+        a.X = np.asarray(a.X.todense(), dtype=np.float32)
+        gc.collect()
+
+    partial_marker = Path(checkpoint_dir, "_partial_epoch") if checkpoint_dir is not None else None
+    is_partial_checkpoint = partial_marker is not None and partial_marker.exists()
+
+    if checkpoint_dir is not None and Path(checkpoint_dir).exists() and not is_partial_checkpoint:
         import torch
 
         map_location = None if torch.cuda.is_available() else torch.device("cpu")
@@ -249,13 +274,31 @@ def build_scooby_embedding_scpoli(
         )
         model = scPoli.load(checkpoint_dir, adata=a, map_location=map_location)
     else:
-        model = scPoli(
-            adata=a,
-            condition_keys=condition_key,
-            cell_type_keys=cell_type_key,
-            recon_loss=recon_loss,
-            latent_dim=latent_dim,
-        )
+        completed_epochs = 0
+        if is_partial_checkpoint:
+            import torch
+
+            completed_epochs = int(partial_marker.read_text())
+            map_location = None if torch.cuda.is_available() else torch.device("cpu")
+            logger.info(
+                "Scooby scPoli embedding: resuming from a PARTIAL checkpoint at %s (%d/%d epochs "
+                "already done), warm-starting from these weights with a fresh optimizer for the "
+                "remaining epochs (not an exact optimizer/scheduler resume like the RNA training "
+                "loop's, since scPoliTrainer has no state_dict for that, but far better than losing "
+                "the whole run)", checkpoint_dir, completed_epochs, n_epochs,
+            )
+            model = scPoli.load(checkpoint_dir, adata=a, map_location=map_location)
+        else:
+            model = scPoli(
+                adata=a,
+                condition_keys=condition_key,
+                cell_type_keys=cell_type_key,
+                recon_loss=recon_loss,
+                latent_dim=latent_dim,
+            )
+
+        remaining_epochs = max(1, n_epochs - completed_epochs)
+        remaining_pretraining_epochs = max(0, pretraining_epochs - completed_epochs)
         early_stopping_kwargs = (
             {
                 "early_stopping_metric": "val_prototype_loss",
@@ -269,16 +312,55 @@ def build_scooby_embedding_scpoli(
             if early_stopping
             else None
         )
-        train_kwargs: dict[str, Any] = {"n_epochs": n_epochs, "pretraining_epochs": pretraining_epochs}
+        train_kwargs: dict[str, Any] = {"n_epochs": remaining_epochs, "pretraining_epochs": remaining_pretraining_epochs}
         if early_stopping_kwargs is not None:
             train_kwargs["early_stopping_kwargs"] = early_stopping_kwargs
-        model.train(**train_kwargs)
+
+        if checkpoint_dir is not None:
+            from scarches.trainers.scpoli.trainer import scPoliTrainer
+
+            _original_on_epoch_end = scPoliTrainer.on_epoch_end
+            _original_on_iteration = scPoliTrainer.on_iteration
+            _checkpoint_every_n_epochs = 5
+            _checkpoint_every_n_iters = 2000
+
+            def _on_epoch_end_with_checkpoint(trainer_self):
+                _original_on_epoch_end(trainer_self)
+                if (trainer_self.epoch + 1) % _checkpoint_every_n_epochs == 0:
+                    real_epoch = completed_epochs + trainer_self.epoch + 1
+                    logger.info(
+                        "Scooby scPoli embedding: periodic checkpoint at epoch %d/%d -> %s",
+                        real_epoch, n_epochs, checkpoint_dir,
+                    )
+                    model.save(checkpoint_dir, overwrite=True)
+                    Path(checkpoint_dir, "_partial_epoch").write_text(str(real_epoch))
+
+            def _on_iteration_with_checkpoint(trainer_self, batch_data):
+                _original_on_iteration(trainer_self, batch_data)
+                if trainer_self.iter > 0 and trainer_self.iter % _checkpoint_every_n_iters == 0:
+                    logger.info(
+                        "Scooby scPoli embedding: mid-epoch checkpoint at epoch %d (iter %d/%d) -> %s "
+                        "(warm-start snapshot, not a completed-epoch boundary)",
+                        completed_epochs + trainer_self.epoch, trainer_self.iter, trainer_self.iters_per_epoch, checkpoint_dir,
+                    )
+                    model.save(checkpoint_dir, overwrite=True)
+                    Path(checkpoint_dir, "_partial_epoch").write_text(str(completed_epochs + trainer_self.epoch))
+
+            scPoliTrainer.on_epoch_end = _on_epoch_end_with_checkpoint
+            scPoliTrainer.on_iteration = _on_iteration_with_checkpoint
+            try:
+                model.train(**train_kwargs)
+            finally:
+                scPoliTrainer.on_epoch_end = _original_on_epoch_end
+                scPoliTrainer.on_iteration = _original_on_iteration
+        else:
+            model.train(**train_kwargs)
 
         if checkpoint_dir is not None:
             logger.info("Training done. saving model to %s.", checkpoint_dir)
             model.save(checkpoint_dir, overwrite=True)
-
-    import scipy.sparse
+            if Path(checkpoint_dir, "_partial_epoch").exists():
+                Path(checkpoint_dir, "_partial_epoch").unlink()
 
     if scipy.sparse.issparse(a.X):
         a.X = np.asarray(a.X.todense())
@@ -427,6 +509,38 @@ def filter_val(df):
     return train_ds, val_ds
 
 
+_PROGRESS_FILENAME = "_cellink_progress.json"
+
+
+def _write_training_progress(checkpoint_dir: str, *, completed_epochs: int, num_epochs: int) -> None:
+    """Persist the true number of fully-completed epochs alongside a checkpoint.
+
+    This is the single source of truth for resume, not the epoch number
+    embedded in the checkpoint directory's own name (``scooby_epoch_{epoch}_
+    {step}_...``): that name-based scheme can drift from the real count
+    across a long chain of resumes, so every checkpoint carries its own
+    verified count directly, immune to directory-name drift, a mid-chain
+    kill/resubmit, or any future bug in the naming scheme.
+    """
+    with open(Path(checkpoint_dir) / _PROGRESS_FILENAME, "w") as f:
+        json.dump({"completed_epochs": completed_epochs, "num_epochs": num_epochs}, f)
+
+
+def _read_training_progress(checkpoint_dir: str) -> int | None:
+    """Read back the true completed-epoch count written by ``_write_training_progress``.
+
+    Returns ``None`` (not 0) when absent, so callers can distinguish an
+    old-style checkpoint that predates this tracking, which should fall back
+    to a best-effort directory-name parse, from a run that has genuinely
+    completed zero epochs.
+    """
+    path = Path(checkpoint_dir) / _PROGRESS_FILENAME
+    if not path.is_file():
+        return None
+    with open(path) as f:
+        return int(json.load(f)["completed_epochs"])
+
+
 def _run_scooby_training_loop(
     *,
     scooby_model,
@@ -445,8 +559,21 @@ def _run_scooby_training_loop(
     mode: str,
     save_every_n_steps: int = 1000,
     log_with: str | None = None,
+    mixed_precision: str | None = "bf16",
+    gradient_accumulation_steps: int = 1,
+    resume_from_checkpoint: str | None = None,
 ):
     """Shared Accelerate-based training loop for both RNA-only and multiome scooby fine-tuning.
+
+    ``resume_from_checkpoint``: path to a directory previously written by
+    ``accelerator.save_state`` (e.g. one of this same function's own
+    ``{output_dir}/scooby_epoch_{epoch}_{i}_{run_name}`` checkpoints). Loaded
+    via ``accelerator.load_state`` after ``accelerator.prepare``, the correct
+    order since model/optimizer/scheduler must already be wrapped before
+    their states can be restored into the wrapped objects. Model weights,
+    optimizer state, scheduler state, and dataloader sampler position are all
+    restored, a genuine continuation rather than a fresh run seeded with old
+    weights.
     """
     import torch
     import torch.nn as nn
@@ -457,10 +584,18 @@ def _run_scooby_training_loop(
     from scooby.utils.utils import add_weight_decay, evaluate, poisson_multinomial_torch
 
     ddp_kwargs = DistributedDataParallelKwargs(static_graph=True)
-    accelerator = Accelerator(log_with=log_with, kwargs_handlers=[ddp_kwargs], step_scheduler_with_optimizer=False)
+    accelerator = Accelerator(
+        log_with=log_with, kwargs_handlers=[ddp_kwargs], step_scheduler_with_optimizer=False,
+        mixed_precision=mixed_precision, gradient_accumulation_steps=gradient_accumulation_steps,
+    )
     device = accelerator.device
 
-    num_steps = (45_000 * num_epochs) // training_loader.batch_size
+    # In real optimizer-update terms (matching warmup_steps' own units, and
+    # the reference's "warmed up over the first 1,000 steps" of real steps,
+    # not raw dataloader batches), divided by gradient_accumulation_steps so
+    # the LR schedule decays over the intended number of real updates rather
+    # than scaling down with the micro-batch count.
+    num_steps = (45_000 * num_epochs) // (training_loader.batch_size * gradient_accumulation_steps)
     parameters = add_weight_decay(scooby_model, lr=lr, weight_decay=wd)
     optimizer = torch.optim.AdamW(parameters)
 
@@ -472,42 +607,85 @@ def _run_scooby_training_loop(
     scooby_model, optimizer, scheduler, training_loader, val_loader = accelerator.prepare(
         scooby_model, optimizer, scheduler, training_loader, val_loader
     )
+
+    starting_epoch = 0
+    if resume_from_checkpoint:
+        accelerator.load_state(resume_from_checkpoint)
+        logger.info("Scooby training: resumed from %s", resume_from_checkpoint)
+        progress = _read_training_progress(resume_from_checkpoint)
+        if progress is not None:
+            starting_epoch = progress
+            logger.info(
+                "Scooby training: resuming at true epoch %d/%d (read from this checkpoint's own "
+                "%s, verified, not inferred)", starting_epoch, num_epochs, _PROGRESS_FILENAME,
+            )
+        else:
+            m = re.search(r"scooby_epoch_(\d+)_\d+_", Path(resume_from_checkpoint).name)
+            if m:
+                starting_epoch = int(m.group(1))
+                logger.warning(
+                    "Scooby training: %s not found in %r (an old-style checkpoint predating real "
+                    "progress-tracking), falling back to parsing epoch %d from the directory name. "
+                    "This inferred value is only as trustworthy as the name itself, which can drift "
+                    "from the true epoch count across chained resumes. Every checkpoint saved from "
+                    "this run onward carries its own verified count and will not need this fallback.",
+                    _PROGRESS_FILENAME, Path(resume_from_checkpoint).name, starting_epoch,
+                )
+            else:
+                logger.warning(
+                    "Scooby training: could not parse an epoch number from checkpoint name %r either, "
+                    "restarting the epoch loop at 0 (real optimizer/scheduler state is still correctly "
+                    "restored regardless; only the redundant-raw-pass bound above is lost for this resume)",
+                    Path(resume_from_checkpoint).name,
+                )
     if log_with:
         accelerator.init_trackers("scooby", init_kwargs={"wandb": {"name": run_name}})
     loss_fn = poisson_multinomial_torch
 
-    for epoch in range(num_epochs):
+    completed_epochs = starting_epoch
+    while completed_epochs < num_epochs:
+        epoch = completed_epochs
         for i, (inputs, rc_augs, targets, cell_emb_idx) in enumerate(tqdm.tqdm(training_loader)):
-            inputs = inputs.permute(0, 2, 1).to(device, non_blocking=True)
-            targets = targets.to(device, non_blocking=True)
-            for rc_aug_idx in rc_augs.nonzero():
-                rc_aug_idx = rc_aug_idx[0]
-                flipped = torch.flip(targets[rc_aug_idx].unsqueeze(0), (1, -3))
-                targets[rc_aug_idx] = fix_rev_comp_fn(flipped)[0]
-            optimizer.zero_grad()
-            with torch.autocast(accelerator.device.type):
-                outputs = scooby_model(inputs, cell_emb_idx)
-                loss = loss_fn(outputs, targets, total_weight=total_weight)
-                if log_with:
-                    accelerator.log({"loss": loss})
-            accelerator.backward(loss)
-            accelerator.clip_grad_norm_(scooby_model.parameters(), clip_global_norm)
-            if log_with:
-                accelerator.log({"learning_rate": scheduler.get_last_lr()[0]})
-            optimizer.step()
-            scheduler.step()
+            with accelerator.accumulate(scooby_model):
+                inputs = inputs.permute(0, 2, 1).to(device, non_blocking=True)
+                targets = targets.to(device, non_blocking=True)
+                for rc_aug_idx in rc_augs.nonzero():
+                    rc_aug_idx = rc_aug_idx[0]
+                    flipped = torch.flip(targets[rc_aug_idx].unsqueeze(0), (1, -3))
+                    targets[rc_aug_idx] = fix_rev_comp_fn(flipped)[0]
+
+                with accelerator.autocast():
+                    outputs = scooby_model(inputs, cell_emb_idx)
+                    loss = loss_fn(outputs, targets, total_weight=total_weight)
+                    if log_with:
+                        accelerator.log({"loss": loss})
+                accelerator.backward(loss)
+                if accelerator.sync_gradients:
+                    accelerator.clip_grad_norm_(scooby_model.parameters(), clip_global_norm)
+
+                optimizer.step()
+                if accelerator.sync_gradients:
+                    scheduler.step()
+                    if log_with:
+                        accelerator.log({"learning_rate": scheduler.get_last_lr()[0]})
+                optimizer.zero_grad()
             if i % eval_every_n == 0:
                 evaluate(accelerator, scooby_model, val_loader, mode=mode, stop_idx=0)
                 scooby_model.train()
             if (i % save_every_n_steps == 0 and epoch != 0) or (i % (2 * save_every_n_steps) == 0 and epoch == 0 and i != 0):
-                accelerator.save_state(output_dir=f"{output_dir}/scooby_epoch_{epoch}_{i}_{run_name}")
-        logger.info("Scooby training: completed epoch %d/%d", epoch + 1, num_epochs)
+                ckpt_dir = f"{output_dir}/scooby_epoch_{epoch}_{i}_{run_name}"
+                accelerator.save_state(output_dir=ckpt_dir)
+
+                _write_training_progress(ckpt_dir, completed_epochs=completed_epochs, num_epochs=num_epochs)
+        completed_epochs += 1
+        logger.info("Scooby training: completed epoch %d/%d (true, persisted count)", completed_epochs, num_epochs)
 
     final_dir = f"{output_dir}/scooby_final_{run_name}"
     accelerator.save_state(output_dir=final_dir)
+    _write_training_progress(final_dir, completed_epochs=completed_epochs, num_epochs=num_epochs)
     if log_with:
         accelerator.end_training()
-    logger.info("Scooby training complete. Final state: %s", final_dir)
+    logger.info("Scooby training complete. Final state: %s (true epoch count: %d/%d)", final_dir, completed_epochs, num_epochs)
     return final_dir
 
 
@@ -544,10 +722,19 @@ def train_scooby(
     num_workers: int = 8,
     log_with: str | None = None,
     runner: ScoobyRunner | None = None,
+    use_transform_borzoi_emb: bool = False,
+    mixed_precision: str | None = "bf16",
+    gradient_accumulation_steps: int = 1,
+    resume_from_checkpoint: str | None = None,
 ) -> str:
     """RNA-only Scooby fine-tuning: cell embedding + DNA sequence -> per-cell
     RNA coverage. No genotype required (see module docstring).
 
+    ``use_transform_borzoi_emb`` defaults to False: this flag is only correct
+    as True when resuming a training run whose backbone has already
+    converged without it (the reference's own `scooby_reproducibility`
+    resume script uses True for exactly that reason). 
+
     Parameters
     ----------
     rna_plus_path, rna_minus_path : str
@@ -594,7 +781,7 @@ def train_scooby(
         n_tracks=num_tracks,
         return_center_bins_only=True,
         disable_cache=True,
-        use_transform_borzoi_emb=True,
+        use_transform_borzoi_emb=use_transform_borzoi_emb,
     )
     scooby_model = get_lora(scooby_model, train=True, lora_config=runner.lora_config)
 
@@ -619,6 +806,8 @@ def train_scooby(
         num_epochs=num_epochs, lr=lr, wd=wd, warmup_steps=warmup_steps, clip_global_norm=clip_global_norm,
         eval_every_n=eval_every_n, total_weight=total_weight, output_dir=output_dir, run_name=run_name,
         fix_rev_comp_fn=fix_rev_comp_rna, mode="rna", log_with=log_with,
+        mixed_precision=mixed_precision, gradient_accumulation_steps=gradient_accumulation_steps,
+        resume_from_checkpoint=resume_from_checkpoint,
     )
 
 
@@ -657,18 +846,23 @@ def train_scooby_multiome(
     num_workers: int = 8,
     log_with: str | None = None,
     runner: ScoobyRunner | None = None,
+    use_transform_borzoi_emb: bool = False,
+    mixed_precision: str | None = "bf16",
+    gradient_accumulation_steps: int = 1,
+    resume_from_checkpoint: str | None = None,
 ) -> str:
     """RNA+ATAC multiome Scooby fine-tuning: cell embedding + DNA sequence ->
-    per-cell RNA coverage AND ATAC accessibility, jointly. 
+    per-cell RNA coverage AND ATAC accessibility, jointly.
 
     Parameters
     ----------
     atac_path : str
         Path to an ATAC h5ad with ``obsm["insertion"]``.
     normalize_atac : bool
-        Scale ATAC coverage (x0.05) for training stability, 
+        Scale ATAC coverage (x0.05) for training stability,
         reference default for multiome training.
-    Other parameters : see :func:`train_scooby`.
+    Other parameters : see :func:`train_scooby` (``use_transform_borzoi_emb``
+        defaults to False here for the same cold-start reason).
 
     Returns
     -------
@@ -701,7 +895,7 @@ def train_scooby_multiome(
         n_tracks=num_tracks,
         return_center_bins_only=True,
         disable_cache=True,
-        use_transform_borzoi_emb=True,
+        use_transform_borzoi_emb=use_transform_borzoi_emb,
     )
     scooby_model = get_lora(scooby_model, train=True, lora_config=runner.lora_config)
 
@@ -726,6 +920,8 @@ def train_scooby_multiome(
         num_epochs=num_epochs, lr=lr, wd=wd, warmup_steps=warmup_steps, clip_global_norm=clip_global_norm,
         eval_every_n=eval_every_n, total_weight=total_weight, output_dir=output_dir, run_name=run_name,
         fix_rev_comp_fn=fix_rev_comp_multiome, mode="multiome", log_with=log_with,
+        mixed_precision=mixed_precision, gradient_accumulation_steps=gradient_accumulation_steps,
+        resume_from_checkpoint=resume_from_checkpoint,
     )
 
 
diff --git a/src/cellink/tl/external/_scprs.py b/src/cellink/tl/external/_scprs.py
index adc62c8..fa53a82 100644
--- a/src/cellink/tl/external/_scprs.py
+++ b/src/cellink/tl/external/_scprs.py
@@ -116,7 +116,7 @@ def _normalise_gwas_df(gwas_df: pd.DataFrame) -> pd.DataFrame:
     has_snp = any(c in df.columns for c in ["SNP", "rsid", "ID", "variant_id"])
     if not has_snp:
         if "CHR" in df.columns and "BP" in df.columns:
-            logger.info("No SNP/rsid column — synthesising variant IDs as 'CHR:BP'.")
+            logger.info("No SNP/rsid column, synthesising variant IDs as 'CHR:BP'.")
             df["SNP"] = (
                 df["CHR"].astype(str).str.replace("chr", "", regex=False)
                 + ":"
@@ -172,7 +172,7 @@ def _annotate_gwas_with_bim_ids(
 
     if gwas_has_rsids and bim_has_rsids:
         logger.info(
-            "SNP ID annotation: GWAS and bim both use rsIDs — skipping "
+            "SNP ID annotation: GWAS and bim both use rsIDs, skipping "
             "CHR:BP lookup (direct rsID matching will be used by PLINK)."
         )
         return gwas_df
@@ -195,8 +195,8 @@ def _annotate_gwas_with_bim_ids(
     )
     if match_pct < 20.0:
         logger.warning(
-            f"Only {match_pct:.1f}% of GWAS variants matched by position — "
-            "this typically means a genome-build mismatch between the GWAS "
+            f"Only {match_pct:.1f}% of GWAS variants matched by position. "
+            "This typically means a genome-build mismatch between the GWAS "
             "(e.g. GRCh38) and the bfile (e.g. GRCh37).  "
             "Fetch the GWAS in the same build as your bfile, or supply a "
             "matching reference panel."
@@ -391,7 +391,7 @@ def _build_knn_graph(
     elif "X_pca" in adata.obsm:
         atac_embed = adata.obsm["X_pca"]
     else:
-        logger.info("No LSI/PCA found — computing PCA.")
+        logger.info("No LSI/PCA found, computing PCA.")
         import scipy.sparse as sp
 
         X_bin = (adata.X > 0).astype(float)
@@ -605,10 +605,10 @@ def write_slurm_array_job(
 
     1. Pre-generates per-cell filtered GWAS files in ``prs_dir`` (Python step,
        requires AnnData, must be run once before cluster submission).
-    2. Writes ``{output_dir}/{job_name}_commands.txt`` — one shell command per
+    2. Writes ``{output_dir}/{job_name}_commands.txt``: one shell command per
        cell, where each command runs PLINK clumping + scoring for all r²
        thresholds for that cell.
-    3. Writes ``{output_dir}/{job_name}_array.sh`` — a SLURM array job script
+    3. Writes ``{output_dir}/{job_name}_array.sh``: a SLURM array job script
        where task ``$SLURM_ARRAY_TASK_ID`` executes the corresponding line of
        the commands file.
     4. Logs the ``sbatch`` command and the follow-up Python call to run after
@@ -656,7 +656,7 @@ def write_slurm_array_job(
         ``slurm_cpus`` so the CPU allocation is fully utilised.
     slurm_partition : str, optional
         ``#SBATCH --partition`` value (``-p`` flag).  If ``None``, the line
-        is omitted — useful on clusters where the default partition is fine or
+        is omitted, useful on clusters where the default partition is fine or
         where ``-q``/``--qos`` selects the queue instead.
     slurm_qos : str, optional
         ``#SBATCH -q`` value.  If ``None``, the line is omitted.
@@ -760,7 +760,7 @@ def write_slurm_array_job(
             "# Read the command for this array task (tasks are 0-indexed, sed is 1-indexed)",
             'CMD=$(sed -n "$((SLURM_ARRAY_TASK_ID + 1))p" "$COMMANDS_FILE")',
             'if [ -z "$CMD" ]; then',
-            '    echo "Task $SLURM_ARRAY_TASK_ID: no command (cell skipped — no peak SNPs)."; exit 0',
+            '    echo "Task $SLURM_ARRAY_TASK_ID: no command (cell skipped, no peak SNPs)."; exit 0',
             "fi",
             "",
             'echo "Task $SLURM_ARRAY_TASK_ID starting: $CMD"',
@@ -786,7 +786,7 @@ def write_slurm_array_job(
         f"  Submit with:     {submit_cmd}\n"
         f"\n"
         f"  Memory guidance:\n"
-        f"    Current default: {slurm_mem} — suitable for 1000G EUR (~500 samples, ~8M variants).\n"
+        f"    Current default: {slurm_mem}, suitable for 1000G EUR (~500 samples, ~8M variants).\n"
         f"    For UK Biobank scale (~400k samples): increase to 64G or more.\n"
         f"\n"
         f"  Timing guidance:\n"
@@ -976,7 +976,7 @@ def run_scprs(
     L_np = nx.normalized_laplacian_matrix(graph).toarray()
 
     if "cuda" in device and not torch.cuda.is_available():
-        logger.warning("CUDA not available — falling back to CPU.")
+        logger.warning("CUDA not available, falling back to CPU.")
         device = "cpu"
     dev = torch.device(device)
 
diff --git a/src/cellink/tl/external/_seismic_torch.py b/src/cellink/tl/external/_seismic_torch.py
index d2a6bd9..03e3212 100644
--- a/src/cellink/tl/external/_seismic_torch.py
+++ b/src/cellink/tl/external/_seismic_torch.py
@@ -239,7 +239,7 @@ def run_seismic_torch(
     ----------
     adata : AnnData
         Single-cell data. ``species="mouse"`` is not auto-translated (unlike
-        the R backend) — pass an adata already indexed by the target species'
+        the R backend), pass an adata already indexed by the target species'
         gene symbols if you need cross-species gene mapping.
     magma_file : str or Path
         MAGMA ``.genes.out`` file with columns ``magma_gene_col``/``magma_z_col``.
@@ -302,8 +302,8 @@ def run_seismic_torch(
     shared_genes = scores_df.index.intersection(magma_df.index)
     if len(shared_genes) < 200:
         raise ValueError(
-            f"Only {len(shared_genes)} genes shared between expression data and MAGMA output "
-            "— check that gene identifiers match (gene symbols vs Ensembl IDs)."
+            f"Only {len(shared_genes)} genes shared between expression data and MAGMA output. "
+            "Check that gene identifiers match (gene symbols vs Ensembl IDs)."
         )
     logger.info(f"run_seismic_torch: {len(shared_genes)} genes shared with MAGMA output")
 
diff --git a/src/cellink/tl/external/_sldsc_utils.py b/src/cellink/tl/external/_sldsc_utils.py
index 5881d73..b788a04 100644
--- a/src/cellink/tl/external/_sldsc_utils.py
+++ b/src/cellink/tl/external/_sldsc_utils.py
@@ -730,7 +730,7 @@ def get_magma_gene_loc(
     See Also
     --------
     generate_gene_coord_file : Headed GENE/CHR/START/END format for S-LDSC.
-    run_magma_annotate : MAGMA Step I — uses the gene location file.
+    run_magma_annotate : MAGMA Step I, uses the gene location file.
     """
     from cellink.resources._utils import get_data_home
 
@@ -762,7 +762,7 @@ def get_magma_gene_loc(
     gene_loc = gene_loc.drop_duplicates(subset=["gene"], keep="first")
     gene_loc = gene_loc.sort_values(["chrom", "start"])
 
-    # Write headless (no column names) — MAGMA format
+    # Write headless (no column names), MAGMA format
     gene_loc.to_csv(out_path, sep="\t", index=False, header=False)
     logger.info("Wrote %d gene locations to %s", len(gene_loc), out_path)
     return str(out_path)
diff --git a/src/cellink/tl/external/_tensorqtl.py b/src/cellink/tl/external/_tensorqtl.py
index 7db535c..4f4f973 100644
--- a/src/cellink/tl/external/_tensorqtl.py
+++ b/src/cellink/tl/external/_tensorqtl.py
@@ -75,6 +75,186 @@ def read_tensorqtl_results(
     return results
 
 
+def build_known_cis_eqtls_from_tensorqtl(
+    tensorqtl_parquet_path: str,
+    gene_names: list[str],
+    max_snps_per_gene: int = 1,
+    pval_threshold: float | None = None,
+) -> pd.DataFrame:
+    """
+    Build a ``known_cis_eqtls`` annotation (variants x genes, binary) from a
+    completed TensorQTL nominal cis-scan, for use as a fine-mapping prior.
+
+    Parameters
+    ----------
+    tensorqtl_parquet_path : str
+        Path to a TensorQTL nominal cis-scan parquet file with ``gene``,
+        ``variant_id``, and ``pval`` columns.
+    gene_names : list[str]
+        Genes to build the annotation for.
+    max_snps_per_gene : int, default=1
+        Number of lowest-p-value variant(s) to select per gene.
+    pval_threshold : float, optional
+        If given, only variants with ``pval <= pval_threshold`` are eligible.
+
+    Returns
+    -------
+    pd.DataFrame
+        Binary variant x gene matrix, 1 where that variant is the selected
+        cis-eQTL for that gene.
+    """
+    df = pd.read_parquet(tensorqtl_parquet_path, columns=["gene", "variant_id", "pval"])
+    df = df[df["gene"].isin(set(gene_names))]
+    if pval_threshold is not None:
+        df = df[df["pval"] <= pval_threshold]
+    df = df.sort_values("pval").groupby("gene", sort=False).head(max_snps_per_gene)
+    if df.empty:
+        raise ValueError(
+            f"No cis-eQTL pairs survived filtering from {tensorqtl_parquet_path} "
+            f"(pval_threshold={pval_threshold}); cannot build known_cis_eqtls."
+        )
+    snps = df["variant_id"].unique()
+    genes_with_hits = df["gene"].unique()
+    known = pd.DataFrame(0, index=snps, columns=genes_with_hits, dtype=int)
+    known.values[
+        pd.Index(snps).get_indexer(df["variant_id"]),
+        pd.Index(genes_with_hits).get_indexer(df["gene"]),
+    ] = 1
+    return known
+
+
+def _map_susie_with_prior_weights(
+    genotype_df: pd.DataFrame,
+    variant_df: pd.DataFrame,
+    phenotype_df: pd.DataFrame,
+    phenotype_pos_df: pd.DataFrame,
+    covariates_df: pd.DataFrame,
+    prior_weights: dict[str, pd.Series],
+    L: int = 10,
+    window: int = 1000000,
+    max_iter: int = 500,
+    maf_threshold: float = 0,
+    scaled_prior_variance: float = 0.2,
+    coverage: float = 0.95,
+    min_abs_corr: float = 0.5,
+    estimate_residual_variance: bool = True,
+    estimate_prior_variance: bool = True,
+    tol: float = 1e-3,
+) -> tuple[pd.DataFrame, dict]:
+    """SuSiE fine-mapping with a per-variant prior weight instead of SuSiE's
+    default uniform prior over variants in the cis-window. Uses the same
+    per-phenotype data preparation as ``tensorqtl.susie.map()`` (monomorphic
+    + MAF filtering, covariate residualization); the one real difference
+    is that ``susie.susie()`` is called with ``prior_weights`` supplied
+    (``tensorqtl.susie.map()`` itself never exposes this, see the call site
+    in ``run_tensorqtl``'s ``cis_susie`` branch for why).
+
+    Parameters
+    ----------
+    prior_weights : dict[str, pd.Series]
+        Keyed by phenotype_id; each value is a Series of raw (unnormalized)
+        per-variant weights indexed by variant_id, covering (at least) that
+        phenotype's own cis-window. Variants in the window missing from the
+        Series are given the Series' own median as a neutral (not zero)
+        weight. A hard zero would make them structurally unselectable
+        regardless of signal, a stronger claim than a missing score
+        justifies. Weights are renormalized to sum to 1 automatically.
+
+    Returns
+    -------
+    Same shape as ``tensorqtl.susie.map(..., summary_only=False)``:
+    ``(susie_summary_df, {phenotype_id: {pip, sets, converged, elbo, niter, lbf_variable}})``.
+    """
+    import torch
+    from tensorqtl import genotypeio
+    from tensorqtl.core import Residualizer, calculate_maf, impute_mean
+    from tensorqtl.susie import susie as susie_fit
+
+    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+    residualizer = Residualizer(torch.tensor(covariates_df.values, dtype=torch.float32).to(device))
+    genotype_ix = np.array([genotype_df.columns.tolist().index(i) for i in phenotype_df.columns])
+    genotype_ix_t = torch.from_numpy(genotype_ix).to(device)
+
+    igc = genotypeio.InputGeneratorCis(genotype_df, variant_df, phenotype_df, phenotype_pos_df, window=window)
+    if igc.n_phenotypes == 0:
+        raise ValueError("No valid phenotypes found.")
+
+    susie_summary = []
+    susie_res = {}
+    for phenotype, genotypes, genotype_range, phenotype_id in igc.generate_data(verbose=True):
+        genotypes_t = torch.tensor(genotypes, dtype=torch.float).to(device)
+        genotypes_t = genotypes_t[:, genotype_ix_t]
+        impute_mean(genotypes_t)
+        variant_ids = variant_df.index[genotype_range[0]:genotype_range[-1] + 1].rename("variant_id")
+
+        mask_t = ~(genotypes_t == genotypes_t[:, [0]]).all(1)
+        if maf_threshold > 0:
+            maf_t = calculate_maf(genotypes_t)
+            mask_t &= maf_t >= maf_threshold
+        if mask_t.any():
+            genotypes_t = genotypes_t[mask_t]
+            mask = mask_t.cpu().numpy().astype(bool)
+            variant_ids = variant_ids[mask]
+            genotype_range = genotype_range[mask]
+        if genotypes_t.shape[0] == 0:
+            logger.warning(f"skipping {phenotype_id} (no valid variants)")
+            continue
+
+        if phenotype_id not in prior_weights:
+            raise KeyError(
+                f"prior_weights has no entry for phenotype_id={phenotype_id!r}; "
+                f"provide a Series for every phenotype passed through cis_output, or omit it from cis_output."
+            )
+        pw = prior_weights[phenotype_id]
+        aligned = pw.reindex(variant_ids)
+        n_scored = aligned.notna().sum()
+        fallback = pw.median()
+        aligned = aligned.fillna(fallback)
+        logger.info(
+            f"{phenotype_id}: {genotypes_t.shape[0]} variants survive monomorphic+MAF filter "
+            f"({genotypes_t.shape[1]} samples); {n_scored} of those have a real prior weight, "
+            f"remaining {len(variant_ids) - n_scored} given the pool median ({fallback:.4g}) as a neutral prior"
+        )
+        prior_weights_t = torch.tensor(aligned.values / aligned.values.sum(), dtype=torch.float32).to(device)
+
+        phenotype_t = torch.tensor(phenotype, dtype=torch.float).to(device)
+        genotypes_res_t = residualizer.transform(genotypes_t)
+        phenotype_res_t = residualizer.transform(phenotype_t.reshape(1, -1))
+
+        res = susie_fit(
+            genotypes_res_t.T, phenotype_res_t.T, L=L, scaled_prior_variance=scaled_prior_variance,
+            prior_weights=prior_weights_t, coverage=coverage, min_abs_corr=min_abs_corr,
+            estimate_residual_variance=estimate_residual_variance, estimate_prior_variance=estimate_prior_variance,
+            tol=tol, max_iter=max_iter,
+        )
+
+        af_t = genotypes_t.sum(1) / (2 * genotypes_t.shape[1])
+        res["pip"] = pd.DataFrame({"pip": res["pip"], "af": af_t.cpu().numpy()}, index=variant_ids)
+        logger.info(f"{phenotype_id}: converged={res['converged']}, niter={res['niter']}, "
+                    f"elbo={res['elbo'][-1] if len(res['elbo']) else None}, "
+                    f"max_pip={res['pip']['pip'].max():.6f}, sets['cs'] is None={res['sets']['cs'] is None}")
+        if res["sets"]["cs"] is not None:
+            if res["converged"]:
+                for c in sorted(res["sets"]["cs"], key=lambda x: int(x.replace("L", ""))):
+                    cs = res["sets"]["cs"][c]
+                    p = res["pip"].iloc[cs].copy().reset_index()
+                    p["cs_id"] = c.replace("L", "")
+                    p.insert(0, "phenotype_id", phenotype_id)
+                    susie_summary.append(p)
+                res["lbf_variable"] = res["lbf_variable"][res["sets"]["cs_index"]]
+        copy_keys = ["pip", "sets", "converged", "elbo", "niter", "lbf_variable"]
+        susie_res[phenotype_id] = {k: res[k] for k in copy_keys}
+
+    susie_summary_df = (
+        pd.concat(susie_summary, axis=0).rename(columns={"snp": "variant_id"}).reset_index(drop=True)
+        if susie_summary else pd.DataFrame()
+    )
+    drop_ids = [k for k in susie_res if susie_res[k]["sets"]["cs"] is None]
+    for k in drop_ids:
+        del susie_res[k]
+    return susie_summary_df, susie_res
+
+
 def _run_tensorqtl_python_api(
     mode: str,
     phenotype_df: pd.DataFrame,
@@ -99,6 +279,7 @@ def _run_tensorqtl_python_api(
     fdr: float,
     qvalue_lambda: float,
     seed: int,
+    prior_weights: dict[str, pd.Series] | None = None,
 ) -> pd.DataFrame | tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame] | tuple[dict, pd.DataFrame]:
     """Run TensorQTL via the Python API directly (no subprocess/file export)."""
     try:
@@ -235,18 +416,25 @@ def _run_tensorqtl_python_api(
         pheno_df_sub = phenotype_df.loc[phenotype_ids]
         pheno_pos_df_sub = phenotype_pos_df.loc[phenotype_ids]
 
-        susie_summary, susie_dict = susie.map(
-            genotype_df,
-            variant_df,
-            pheno_df_sub,
-            pheno_pos_df_sub,
-            covariates_df,
-            L=max_effects,
-            window=window,
-            summary_only=False,
-            max_iter=500,
-            maf_threshold=maf_threshold,
-        )
+        if prior_weights is None:
+            susie_summary, susie_dict = susie.map(
+                genotype_df,
+                variant_df,
+                pheno_df_sub,
+                pheno_pos_df_sub,
+                covariates_df,
+                L=max_effects,
+                window=window,
+                summary_only=False,
+                max_iter=500,
+                maf_threshold=maf_threshold,
+            )
+        else:
+            susie_summary, susie_dict = _map_susie_with_prior_weights(
+                genotype_df, variant_df, pheno_df_sub, pheno_pos_df_sub, covariates_df,
+                prior_weights=prior_weights, L=max_effects, window=window, max_iter=500,
+                maf_threshold=maf_threshold,
+            )
         results = (susie_dict, susie_summary)
 
     elif mode == "trans_susie":
@@ -320,6 +508,7 @@ def run_tensorqtl(
     overwrite_covariates_export: bool = True,
     overwrite_phenotype_export: bool = True,
     overwrite_plink_export: bool = True,
+    prior_weights: dict[str, pd.Series] | None = None,
 ) -> pd.DataFrame | tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame] | tuple[dict, pd.DataFrame] | str:
     """
     Run cis- or trans-QTL mapping using TensorQTL on donor-level aggregated expression and genotype data.
@@ -425,6 +614,16 @@ def run_tensorqtl(
         If True, overwrites the phenotype export.
     overwrite_plink_export : bool, default=True
         If True, overwrites the plink export.
+    prior_weights : dict[str, pd.Series], optional
+        ``mode="cis_susie"`` only, requires ``use_python_api=True``. Keyed
+        by phenotype_id; each value is a Series of raw (unnormalized)
+        per-variant prior weights indexed by variant_id, used in place of
+        SuSiE's default flat/uniform prior over variants in the cis-window, 
+        e.g. a sequence-model-predicted regulatory effect size, letting
+        an orthogonal source of information sharpen fine-mapping
+        resolution on credible sets the flat prior alone leaves ambiguous.
+        Variants in a phenotype's window missing from its Series are given
+        that Series' own median as a neutral (not zero) prior. 
 
     Returns
     -------
@@ -445,6 +644,15 @@ def run_tensorqtl(
     if plink_export_kwargs is None:
         plink_export_kwargs = {}
 
+    if prior_weights is not None:
+        if mode != "cis_susie":
+            raise ValueError(f"prior_weights is only supported for mode='cis_susie', got mode={mode!r}.")
+        if not use_python_api:
+            raise ValueError(
+                "prior_weights requires use_python_api=True; it is not expressible through tensorqtl's "
+                "own CLI/subprocess interface, only via cellink's Python-API re-implementation of cis_susie."
+            )
+
     if "X_pca" not in dd.C.obsm:
         logger.info("Calculating PCA.")
         sc.pp.pca(dd.C, n_comps=n_pcs)
@@ -530,6 +738,7 @@ def run_tensorqtl(
             fdr=fdr,
             qvalue_lambda=qvalue_lambda,
             seed=seed,
+            prior_weights=prior_weights,
         )
 
     if run:
diff --git a/src/cellink/tl/external/config/livi.yaml b/src/cellink/tl/external/config/livi.yaml
index fa6c2e0..0961d6a 100644
--- a/src/cellink/tl/external/config/livi.yaml
+++ b/src/cellink/tl/external/config/livi.yaml
@@ -4,8 +4,8 @@
 #   cl.tl.external.configure_livi_runner(livi_root=..., config_path="path/to/livi.yaml")
 #
 # execution_mode options:
-#   python_api  — import LIVI directly (recommended; returns Python objects)
-#   subprocess  — run LIVI's Hydra CLI as a child process (HPC job submission)
+#   python_api  : import LIVI directly (recommended; returns Python objects)
+#   subprocess  : run LIVI's Hydra CLI as a child process (HPC job submission)
 
 execution_mode: python_api
 
diff --git a/tests/test_donor_data.py b/tests/test_donor_data.py
index f70e9c6..56874ca 100644
--- a/tests/test_donor_data.py
+++ b/tests/test_donor_data.py
@@ -6,7 +6,10 @@
 from cellink._core.data_fields import CAnn, DAnn
 from cellink._core.donordata import DonorData
 
-md.set_options(pull_on_update=False)
+if hasattr(md, "set_options"):
+    md.set_options(pull_on_update=False)  # mudata <0.4
+else:
+    md.settings.pull_on_update = False  # mudata >=0.4 (Python >=3.12): set_options -> settings
 
 DATA = Path("tests/data")