fixes

CHANGELOG.md

@@ -12,6 +12,19 @@ and this project adheres to [Semantic Versioning][].
 
 ### Added
 
+- Fixed `SCLINKER_ENHANCER_LINKS_GENOME_BUILD`: the Broad sc-linker Roadmap/ABC
+  enhancer-gene links are GRCh37, not GRCh38. Declaring GRCh38 refused correct
+  GRCh37 setups and, worse, let a GRCh38 panel pass the build check while
+  intersecting hg19 enhancers against hg38 SNPs. Verified from the files: ABC
+  `TargetGeneTSS` matches hg19 exactly for BACH2/CTLA4/FOXP3, and the Roadmap
+  file's largest chr1 coordinate (249,240,000) exceeds GRCh38 chr1's length.
+- `resources.get_eqtl_catalog_credible_sets` and `resources.get_eqtl_catalog_lbf`,
+  exposing the eQTL Catalogue's SuSiE fine-mapping output (per-variant PIPs, and the
+  per-variant log Bayes factors that `tl.coloc_susie` needs for the QTL side -- there
+  was previously no way to obtain these through cellink)
+- `region=` on `resources.get_eqtl_catalog_dataset_associations`, performing a remote
+  tabix range query instead of downloading a whole dataset (a single dataset's
+  summary statistics are ~1.4 GB)
 - Basic tool, preprocessing and plotting functions
 - LIVI donor-level representation learning, sc-linker gene programs, scPRS, gsMap and
   MAGMA wrappers under `cellink.tl.external`, now documented in the API reference
@@ -28,6 +41,12 @@ and this project adheres to [Semantic Versioning][].
 
 ### Fixed
 
+- Both eQTL Catalogue accessors were non-functional: `resources.get_eqtl_catalog_datasets`
+  and `resources.get_eqtl_catalog_dataset_associations` targeted the retired REST API at
+  `https://www.ebi.ac.uk/eqtl/api/v3`, which returns HTTP 410 for every endpoint and
+  version. Both are rewritten against the current FTP/tabix distribution
+  (https://www.ebi.ac.uk/eqtl/Data_access/). `max_pages` is still accepted but ignored
+  with a warning, since the dataset index is no longer paginated
 - `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

docs/api/resources.md

@@ -11,8 +11,10 @@
    resources.get_1000genomes_grch38
    resources.get_dummy_onek1k
    resources.get_onek1k
-   resources.get_eqtl_catalog_dataset_associations
    resources.get_eqtl_catalog_datasets
+   resources.get_eqtl_catalog_dataset_associations
+   resources.get_eqtl_catalog_credible_sets
+   resources.get_eqtl_catalog_lbf
    resources.get_gwas_catalog_studies
    resources.get_gwas_catalog_study
    resources.get_gwas_catalog_study_summary_stats

src/cellink/_core/donordata.py

@@ -194,7 +194,10 @@ def _write_dd(self, f: h5py.File, zarr_path: str | None = None, x_chunks=None):
         f.attrs["var_dims_to_sync"] = self._var_dims_to_sync
 
         for key, value in self.uns.items():
-            f.create_dataset(f"uns/{key}", data=value)
+            try:
+                write_elem(f, f"uns/{key}", value)
+            except (TypeError, NotImplementedError):
+                f.create_dataset(f"uns/{key}", data=value)
 
     def write_h5_dd(self, path: str) -> None:
         """Write the DonorData object to the specified file path.

src/cellink/_core/schema.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import pandas as pd
 import pandera.pandas as pa
 from pandera.pandas import Column, DataFrameSchema
 
@@ -58,11 +59,14 @@ def validate(dd, check_var: bool = True, check_donor_alignment: bool = True) ->
             raise DonorDataSchemaError(f"dd.G.var failed schema validation: {e}") from e
 
     if check_donor_alignment:
-        g_donors = _donor_ids(dd.G, dd.donor_id)
-        c_donors = _donor_ids(dd.C, dd.donor_id)
-        if len(g_donors) != len(c_donors):
+        g_donors = pd.unique(pd.Series(_donor_ids(dd.G, dd.donor_id)))
+        c_donors = pd.unique(pd.Series(_donor_ids(dd.C, dd.donor_id)))
+        if set(g_donors) != set(c_donors):
+            only_g = sorted(set(g_donors) - set(c_donors))
+            only_c = sorted(set(c_donors) - set(g_donors))
             raise DonorDataSchemaError(
-                f"dd.G and dd.C have different donor counts ({len(g_donors)} vs. {len(c_donors)}); "
+                f"dd.G and dd.C cover different donors ({len(g_donors)} vs. {len(c_donors)} "
+                f"unique); only in G: {only_g[:8]}, only in C: {only_c[:8]}. "
                 "DonorData's own construction should never allow this."
             )
         if list(g_donors) != list(c_donors):
@@ -76,10 +80,8 @@ def validate(dd, check_var: bool = True, check_donor_alignment: bool = True) ->
 
 
 def _donor_ids(modality, donor_id: str):
-    from mudata import MuData
-
-    if isinstance(modality, MuData):
-        return modality.obs_names
+    """Donor identifiers for one side of a DonorData, one entry per row of that side.
+    """
     if donor_id in modality.obs.columns:
         return modality.obs[donor_id].to_numpy()
     return modality.obs_names

src/cellink/io/_readwrite.py

@@ -55,7 +55,7 @@ def _read_mudata(group: StorageType, backed: bool = True) -> MuData:
             mods = ModDict()
             gmods = group[k]
             for m in gmods.keys():
-                ad = _read_h5mu_mod(gmods[m], None, True)
+                ad = _read_h5mu_mod(gmods[m], None, False)
                 mods[m] = ad
 
             mod_order = None
@@ -71,7 +71,11 @@ def _read_mudata(group: StorageType, backed: bool = True) -> MuData:
     if "axis" in group.attrs:
         d["axis"] = group.attrs["axis"]
 
-    mu = MuData._init_from_dict_(**d)
+    if hasattr(MuData, "_init_from_dict_"):
+        mu = MuData._init_from_dict_(**d)          # mudata < 0.4
+    else:
+        d["data"] = d.pop("mod", ModDict())        # mudata >= 0.4
+        mu = MuData(**d)
     return mu
 
 
@@ -146,7 +150,11 @@ def _read_anndata(group):
     uns_group = f.get("uns")
     if uns_group:
         for key in uns_group:
-            uns[key] = uns_group[key][()]
+            node = uns_group[key]
+            try:
+                uns[key] = read_elem(node)
+            except Exception:
+                uns[key] = node[()] if hasattr(node, "shape") else node
 
     dd = DonorData(G=G, C=C, donor_id=donor_id, var_dims_to_sync=var_dims_to_sync, uns=uns)
 

src/cellink/resources/__init__.py

@@ -1,7 +1,9 @@
 from ._datasets import get_1000genomes, get_1000genomes_grch38, get_dummy_onek1k, get_onek1k
 from ._gwas_prs_qtl import (
+    get_eqtl_catalog_credible_sets,
     get_eqtl_catalog_dataset_associations,
     get_eqtl_catalog_datasets,
+    get_eqtl_catalog_lbf,
     get_gwas_catalog_studies,
     get_gwas_catalog_study,
     get_gwas_catalog_study_summary_stats,