From 2cc416bb6838d63c08df4a653215c9129aed6cd7 Mon Sep 17 00:00:00 2001 From: Lucas Arnoldt Date: Sun, 6 Sep 2026 03:43:01 +0200 Subject: [PATCH 1/2] fixes --- CHANGELOG.md | 19 ++ docs/api/resources.md | 4 +- src/cellink/_core/donordata.py | 5 +- src/cellink/_core/schema.py | 18 +- src/cellink/io/_readwrite.py | 14 +- src/cellink/resources/__init__.py | 2 + src/cellink/resources/_gwas_prs_qtl.py | 378 ++++++++++++++++++--- src/cellink/tl/external/_scdrs.py | 11 +- src/cellink/tl/external/_sclinker_utils.py | 2 +- 9 files changed, 392 insertions(+), 61 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bf3eb30..95c3946 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/api/resources.md b/docs/api/resources.md index 30e881c..c8f0bc6 100644 --- a/docs/api/resources.md +++ b/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 diff --git a/src/cellink/_core/donordata.py b/src/cellink/_core/donordata.py index 06742f4..9000be4 100644 --- a/src/cellink/_core/donordata.py +++ b/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. diff --git a/src/cellink/_core/schema.py b/src/cellink/_core/schema.py index 64e20ae..78d2973 100644 --- a/src/cellink/_core/schema.py +++ b/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 diff --git a/src/cellink/io/_readwrite.py b/src/cellink/io/_readwrite.py index e63b100..058a95f 100644 --- a/src/cellink/io/_readwrite.py +++ b/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) diff --git a/src/cellink/resources/__init__.py b/src/cellink/resources/__init__.py index e47a7fc..091c7bd 100644 --- a/src/cellink/resources/__init__.py +++ b/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, diff --git a/src/cellink/resources/_gwas_prs_qtl.py b/src/cellink/resources/_gwas_prs_qtl.py index ee2cf39..c277d0b 100644 --- a/src/cellink/resources/_gwas_prs_qtl.py +++ b/src/cellink/resources/_gwas_prs_qtl.py @@ -1,5 +1,9 @@ +import io import logging import re +import shutil +import subprocess +from collections.abc import Sequence from pathlib import Path from typing import Any from urllib.error import URLError @@ -8,7 +12,7 @@ import pandas as pd import requests -from cellink.resources._utils import _cache_df, _to_dataframe, get_data_home +from cellink.resources._utils import _cache_df, _download_file, _to_dataframe, get_data_home def _normalize_build(genome_build: str) -> str: @@ -45,7 +49,20 @@ def _find_candidate_files(html: str) -> list[str]: GWAS_API_BASE = "https://www.ebi.ac.uk/gwas/rest/api/v2" PGS_API_BASE = "https://www.pgscatalog.org/rest" -EQTL_API_BASE = "https://www.ebi.ac.uk/eqtl/api/v3" +EQTL_FTP_BASE = "https://ftp.ebi.ac.uk/pub/databases/spot/eQTL" +EQTL_TABIX_PATHS_URL = ( + "https://raw.githubusercontent.com/eQTL-Catalogue/eQTL-Catalogue-resources/" + "master/tabix/tabix_ftp_paths.tsv" +) +EQTL_SUMSTATS_COLUMNS = ( + "molecular_trait_id", "chromosome", "position", "ref", "alt", "variant", "ma_samples", + "maf", "pvalue", "beta", "se", "type", "ac", "an", "r2", "molecular_trait_object_id", + "gene_id", "median_tpm", "rsid", +) +EQTL_CREDIBLE_SET_COLUMNS = ( + "molecular_trait_id", "gene_id", "cs_id", "variant", "rsid", "cs_size", "pip", + "pvalue", "beta", "se", "z", "cs_min_r2", "region", +) def _fetch( @@ -71,28 +88,37 @@ def _fetch( If the endpoint supports pagination, returns a list of results aggregated across pages. Otherwise, returns the raw JSON response as a dictionary. """ - results = [] + results: list = [] page = 0 - while url: - logging.info(f"Fetching {url}") - r = requests.get(url, params=params) + next_url: str | None = url + next_params = params + + while next_url: + logging.info(f"Fetching {next_url}") + r = requests.get(next_url, params=next_params) r.raise_for_status() data = r.json() + page_items: list | None = None if "_embedded" in data: - for v in data["_embedded"].values(): - results.extend(v) + page_items = [item for v in data["_embedded"].values() for item in v] elif "results" in data: - results.extend(data["results"]) - else: + page_items = data["results"] + + if page_items is None: + if results: + logging.debug(f"No collection payload on page {page}; returning {len(results)} collected items.") + break return data - if "_links" in data: - url = data.get("_links", {}).get("next", {}).get("href") if paginate else None - elif "next" in data: - url = data["next"] if paginate else None - else: - url = None + results.extend(page_items) + + if not paginate: + break + + links = data.get("_links") or {} + next_url = (links.get("next") or {}).get("href") or None if links else data.get("next") or None + next_params = None page += 1 if max_pages and page >= max_pages: @@ -608,35 +634,154 @@ def get_pgs_catalog_score_file( return df +def _eqtl_https(url: str) -> str: + """Rewrite an eQTL Catalogue ``ftp://`` path to its ``https://`` equivalent. + + The catalogue's metadata table publishes ``ftp://ftp.ebi.ac.uk/...`` URLs, but the + same tree is served over HTTPS with byte-range support, which is what htslib needs + for remote tabix queries (and what works from behind proxies that block FTP). + """ + if url.startswith("ftp://"): + return "https://" + url[len("ftp://") :] + return url + + +def _eqtl_tabix_query(url: str, regions: Sequence[str], columns: Sequence[str]) -> pd.DataFrame: + """Run a remote tabix range query against a bgzipped eQTL Catalogue file. + + A single dataset's ``.all.tsv.gz`` is ~1.4 GB, so region-restricted access is the + only practical way to read it; the catalogue ships a ``.tbi`` alongside every + sumstats file for exactly this purpose. + + Parameters + ---------- + url : str + HTTPS URL of the bgzipped, tabix-indexed file. + regions : sequence of str + Regions in tabix syntax, e.g. ``["6:89900000-90300000"]``. Note the catalogue + indexes chromosomes **without** a ``chr`` prefix. + columns : sequence of str + Column names to apply; the catalogue's files carry a header line that is *not* + marked as a comment, so it is absent from the tabix index and cannot be + recovered with ``tabix -H``. + + Returns + ------- + pd.DataFrame + """ + if shutil.which("tabix") is None: + raise RuntimeError( + "The `tabix` executable is required for region-restricted eQTL Catalogue queries " + "but was not found on $PATH. Install htslib (e.g. `conda install -c bioconda htslib`), " + "or call this function without `region=` to download the whole dataset instead." + ) + + frames: list[pd.DataFrame] = [] + for region in regions: + cmd = ["tabix", url, region] + logging.info(f"tabix {url} {region}") + proc = subprocess.run(cmd, capture_output=True, text=True) + if proc.returncode != 0: + raise RuntimeError(f"tabix failed for region {region!r}: {proc.stderr.strip()[:500]}") + if not proc.stdout.strip(): + logging.warning(f"No eQTL Catalogue records returned for region {region!r}.") + continue + frames.append(pd.read_csv(io.StringIO(proc.stdout), sep="\t", header=None, names=list(columns))) + + if not frames: + return pd.DataFrame(columns=list(columns)) + return pd.concat(frames, ignore_index=True) + + def get_eqtl_catalog_datasets( - data_home: str | Path | None = None, max_pages: int | None = None, refresh: bool = False, **params: Any + data_home: str | Path | None = None, + max_pages: int | None = None, + refresh: bool = False, + **params: Any, ) -> pd.DataFrame: """ - Retrieve eQTL catalog datasets and cache locally. + Retrieve the eQTL Catalogue dataset index and cache locally. + + Returns one row per dataset (a study x sample-group x quantification-method + combination), including the FTP paths of its summary statistics, SuSiE credible + sets and log-Bayes-factor files. Parameters ---------- data_home : str or Path, optional Directory to store cached files. Defaults to user data directory. max_pages : int, optional - Maximum number of API pages to fetch. + Ignored. Retained only so that code written against the retired REST API keeps + working; the index is now a single file with no pagination. refresh : bool, default=False - If True, ignore cached data and fetch fresh data. + If True, ignore cached data and re-download the index. **params - Additional query parameters to filter datasets. + Row filters applied to the returned table, as ``column=value`` (or + ``column=[v1, v2]`` for a membership test). Useful columns are + ``study_label``, ``tissue_label``, ``condition_label``, ``quant_method``, + ``dataset_id`` and ``study_id``. String comparisons are case-insensitive. Returns ------- pd.DataFrame - DataFrame containing eQTL catalog datasets metadata. + Columns: ``study_id``, ``dataset_id``, ``study_label``, ``sample_group``, + ``tissue_id``, ``tissue_label``, ``condition_label``, ``sample_size``, + ``quant_method``, ``ftp_path``, ``ftp_cs_path``, ``ftp_lbf_path``. + + Examples + -------- + >>> import cellink as cl + >>> tregs = cl.resources.get_eqtl_catalog_datasets(quant_method="ge", tissue_label="Treg memory") + >>> tregs[["dataset_id", "study_label", "sample_size"]] """ + if max_pages is not None: + logging.warning( + "`max_pages` is ignored: the eQTL Catalogue REST API was retired and the dataset " + "index is now a single unpaginated file." + ) + data_home = get_data_home(data_home) - return _cache_df( - data_home, - "eqtl_datasets.parquet", - refresh, - lambda: _to_dataframe(_fetch(f"{EQTL_API_BASE}/datasets", params=params, max_pages=max_pages)), - ) + + def _fetch_index() -> pd.DataFrame: + logging.info(f"Fetching eQTL Catalogue dataset index from {EQTL_TABIX_PATHS_URL}") + return pd.read_csv(EQTL_TABIX_PATHS_URL, sep="\t") + + df = _cache_df(data_home, "eqtl_datasets.parquet", refresh, _fetch_index) + + for key, value in params.items(): + if key not in df.columns: + raise KeyError(f"'{key}' is not a column of the eQTL Catalogue index. Available: {list(df.columns)}") + col = df[key] + if isinstance(value, (list, tuple, set)): + wanted = {str(v).casefold() for v in value} + df = df[col.astype(str).str.casefold().isin(wanted)] + else: + df = df[col.astype(str).str.casefold() == str(value).casefold()] + + return df.reset_index(drop=True) + + +def _resolve_eqtl_dataset_path( + dataset_id: str, + path_column: str, + data_home: str | Path | None, + refresh: bool, +) -> str: + """Look up one dataset's file URL in the catalogue index.""" + index = get_eqtl_catalog_datasets(data_home=data_home, refresh=refresh) + hit = index[index["dataset_id"] == dataset_id] + if hit.empty: + raise KeyError( + f"dataset_id '{dataset_id}' is not in the eQTL Catalogue index. " + "List available datasets with `cellink.resources.get_eqtl_catalog_datasets()`." + ) + url = hit.iloc[0][path_column] + if not isinstance(url, str) or not url or url.upper() == "NA": + raise ValueError( + f"dataset '{dataset_id}' has no '{path_column}' entry " + "(not every dataset has SuSiE fine-mapping results)." + ) + return _eqtl_https(url) def get_eqtl_catalog_dataset_associations( @@ -644,43 +789,185 @@ def get_eqtl_catalog_dataset_associations( data_home: str | Path | None = None, refresh: bool = False, return_path: bool = False, + region: str | Sequence[str] | None = None, **params: Any, ) -> pd.DataFrame | Path: """ - Retrieve associations for a specific eQTL catalog dataset and cache locally. + Retrieve cis-QTL summary statistics for one eQTL Catalogue dataset. Parameters ---------- dataset_id : str - eQTL catalog dataset ID (e.g., "QTD000319"). + eQTL Catalogue dataset ID (e.g., ``"QTD000625"``, OneK1K Treg memory). data_home : str or Path, optional Directory to store cached files. Defaults to user data directory. refresh : bool, default=False - If True, ignore cached data and fetch fresh data. + If True, ignore cached data and re-fetch. return_path : bool, default=False - If True, return the local cached file path instead of reading it into a DataFrame. + If True, return the local cached file path instead of a DataFrame. Only + meaningful for a whole-dataset download (``region=None``). + region : str or sequence of str, optional + One or more regions in tabix syntax (``"6:89900000-90300000"``), fetched with a + remote range query instead of downloading the file. **Chromosomes are named + without a ``chr`` prefix** and coordinates are GRCh38. Strongly recommended: + a single dataset's full summary statistics are ~1.4 GB. **params - Additional query parameters to pass to the API. + Post-hoc row filters applied to the result as ``column=value``, e.g. + ``gene_id="ENSG00000112182"`` (BACH2) or ``rsid="rs72928038"``. Returns ------- pd.DataFrame or Path - DataFrame of eQTL associations, or Path to the cached parquet file if `return_path=True`. + Summary statistics, or the cached file path when ``return_path=True``. + + Notes + ----- + The eQTL Catalogue REST API this function previously used + (``https://www.ebi.ac.uk/eqtl/api/v3``) was retired and now returns HTTP 410 for + every endpoint and version. Access is via the FTP/tabix distribution described at + https://www.ebi.ac.uk/eqtl/Data_access/. + + Examples + -------- + >>> import cellink as cl + >>> # cis-eQTLs at the BACH2 locus in OneK1K memory Tregs + >>> df = cl.resources.get_eqtl_catalog_dataset_associations( + ... "QTD000625", region="6:89900000-90300000" + ... ) + >>> df.nsmallest(5, "pvalue")[["rsid", "gene_id", "pvalue", "beta"]] """ data_home = get_data_home(data_home) - dest = data_home / f"{dataset_id}_eqtl_associations.parquet" + url = _resolve_eqtl_dataset_path(dataset_id, "ftp_path", data_home, refresh) + + if region is not None: + regions = [region] if isinstance(region, str) else list(region) + df = _eqtl_tabix_query(url, regions, EQTL_SUMSTATS_COLUMNS) + if return_path: + dest = data_home / f"{dataset_id}_eqtl_associations_{'_'.join(regions).replace(':', '-')}.parquet" + df.to_parquet(dest, index=False) + return dest + else: + dest = data_home / f"{dataset_id}.all.tsv.gz" + if not dest.exists() or refresh: + logging.warning( + f"Downloading the complete summary statistics for {dataset_id} (~1 GB or more). " + "Pass `region=` for a remote tabix range query instead." + ) + _download_file(url, dest) + if return_path: + return dest + df = pd.read_csv(dest, sep="\t") + + for key, value in params.items(): + if key not in df.columns: + raise KeyError(f"'{key}' is not a column of the summary statistics. Available: {list(df.columns)}") + if isinstance(value, (list, tuple, set)): + df = df[df[key].isin(list(value))] + else: + df = df[df[key] == value] + + return df.reset_index(drop=True) + - if dest.exists() and not refresh: - return dest if return_path else pd.read_parquet(dest) +def get_eqtl_catalog_credible_sets( + dataset_id: str, + data_home: str | Path | None = None, + refresh: bool = False, + return_path: bool = False, + **params: Any, +) -> pd.DataFrame | Path: + """ + Retrieve SuSiE fine-mapped credible sets for one eQTL Catalogue dataset. - data = _fetch(f"{EQTL_API_BASE}/datasets/{dataset_id}/associations", params=params) - df = _to_dataframe(data) - df.to_parquet(dest) + These are the per-variant posterior inclusion probabilities used to define causal + variants (e.g. the ``PIP >= 0.9`` threshold adopted for variant-effect benchmarks in + the Borzoi and scooby papers), and are the natural input to a colocalization + analysis against a GWAS. + Parameters + ---------- + dataset_id : str + eQTL Catalogue dataset ID (e.g., ``"QTD000625"``). + data_home : str or Path, optional + Directory to store cached files. Defaults to user data directory. + refresh : bool, default=False + If True, ignore cached data and re-download. + return_path : bool, default=False + If True, return the local cached file path instead of a DataFrame. + **params + Row filters applied to the result as ``column=value``, e.g. ``gene_id=...``. + + Returns + ------- + pd.DataFrame or Path + Columns: ``molecular_trait_id``, ``gene_id``, ``cs_id``, ``variant``, ``rsid``, + ``cs_size``, ``pip``, ``pvalue``, ``beta``, ``se``, ``z``, ``cs_min_r2``, + ``region``. + + Examples + -------- + >>> import cellink as cl + >>> cs = cl.resources.get_eqtl_catalog_credible_sets("QTD000625") + >>> cs[cs.pip >= 0.9].head() + """ + data_home = get_data_home(data_home) + url = _resolve_eqtl_dataset_path(dataset_id, "ftp_cs_path", data_home, refresh) + dest = data_home / f"{dataset_id}.credible_sets.tsv.gz" + if not dest.exists() or refresh: + _download_file(url, dest) if return_path: return dest - return df + df = pd.read_csv(dest, sep="\t") + for key, value in params.items(): + if key not in df.columns: + raise KeyError(f"'{key}' is not a column of the credible sets. Available: {list(df.columns)}") + if isinstance(value, (list, tuple, set)): + df = df[df[key].isin(list(value))] + else: + df = df[df[key] == value] + return df.reset_index(drop=True) + + +def get_eqtl_catalog_lbf( + dataset_id: str, + data_home: str | Path | None = None, + refresh: bool = False, + return_path: bool = False, +) -> pd.DataFrame | Path: + """ + Retrieve per-variant SuSiE log Bayes factors for one eQTL Catalogue dataset. + + This is the input :func:`cellink.tl.coloc_susie` expects for the QTL side of a + SuSiE-based colocalization (one column of log Bayes factors per SuSiE component + ``L1..L10``), as opposed to the single-causal-variant approximation used by + :func:`cellink.tl.coloc_abf`. + + Parameters + ---------- + dataset_id : str + eQTL Catalogue dataset ID (e.g., ``"QTD000625"``). + data_home : str or Path, optional + Directory to store cached files. Defaults to user data directory. + refresh : bool, default=False + If True, ignore cached data and re-download. + return_path : bool, default=False + If True, return the local cached file path instead of a DataFrame. Recommended: + these files are ~100 MB compressed and cover every fine-mapped region in the + dataset, so reading a whole one into memory is usually not what you want. + + Returns + ------- + pd.DataFrame or Path + """ + data_home = get_data_home(data_home) + url = _resolve_eqtl_dataset_path(dataset_id, "ftp_lbf_path", data_home, refresh) + dest = data_home / f"{dataset_id}.lbf_variable.txt.gz" + if not dest.exists() or refresh: + _download_file(url, dest) + if return_path: + return dest + return pd.read_csv(dest, sep="\t") if __name__ == "__main__": @@ -708,7 +995,12 @@ def get_eqtl_catalog_dataset_associations( pgs_score_file = get_pgs_catalog_score_file("PGS000043") - eqtl_datasets = get_eqtl_catalog_datasets() + eqtl_datasets = get_eqtl_catalog_datasets(quant_method="ge") print(eqtl_datasets.head()) - eqtl_dataset = get_eqtl_catalog_dataset_associations("QTD000319") + # Region-restricted query: cis-eQTLs at the BACH2 locus in OneK1K memory Tregs. + eqtl_dataset = get_eqtl_catalog_dataset_associations("QTD000625", region="6:89900000-90300000") + print(eqtl_dataset.nsmallest(5, "pvalue")) + + credible_sets = get_eqtl_catalog_credible_sets("QTD000625") + print(credible_sets[credible_sets.pip >= 0.9].head()) diff --git a/src/cellink/tl/external/_scdrs.py b/src/cellink/tl/external/_scdrs.py index 39912c3..8189338 100644 --- a/src/cellink/tl/external/_scdrs.py +++ b/src/cellink/tl/external/_scdrs.py @@ -167,10 +167,13 @@ def run_scdrs( if "X_pca" not in adata.obsm_keys(): logger.info(f"Computing PCA with {n_pcs} components") - sc.pp.highly_variable_genes(adata, n_top_genes=2000) - adata = adata[:, adata.var.highly_variable] - sc.pp.scale(adata, max_value=10) - sc.tl.pca(adata, n_comps=n_pcs) + pca_view = adata.copy() + sc.pp.highly_variable_genes(pca_view, n_top_genes=2000) + pca_view = pca_view[:, pca_view.var.highly_variable].copy() + sc.pp.scale(pca_view, max_value=10) + sc.tl.pca(pca_view, n_comps=n_pcs) + adata.obsm["X_pca"] = pca_view.obsm["X_pca"] + del pca_view adata.X = adata.layers["counts"] diff --git a/src/cellink/tl/external/_sclinker_utils.py b/src/cellink/tl/external/_sclinker_utils.py index 671fccc..0f29699 100644 --- a/src/cellink/tl/external/_sclinker_utils.py +++ b/src/cellink/tl/external/_sclinker_utils.py @@ -1024,7 +1024,7 @@ def genescores_to_100kb_bedgraph( return bedgraphs -SCLINKER_ENHANCER_LINKS_GENOME_BUILD = "GRCh38" +SCLINKER_ENHANCER_LINKS_GENOME_BUILD = "GRCh37" def _check_build_consistency( From 9cf3bf56fa322b44b5e6f160f6a23925f7af0008 Mon Sep 17 00:00:00 2001 From: Lucas Arnoldt Date: Sun, 6 Sep 2026 14:27:31 +0200 Subject: [PATCH 2/2] fixes --- src/cellink/tl/external/_sclinker_utils.py | 11 ++++++----- tests/test_sclinker_build_consistency.py | 5 +++-- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/cellink/tl/external/_sclinker_utils.py b/src/cellink/tl/external/_sclinker_utils.py index 0f29699..fd385b7 100644 --- a/src/cellink/tl/external/_sclinker_utils.py +++ b/src/cellink/tl/external/_sclinker_utils.py @@ -540,11 +540,12 @@ def _query_biomart_and_write_gene_coords(data_dir: Path, genome_build: str = "GR genome_build ``"GRCh37"`` (default) or ``"GRCh38"``. Selects the matching Ensembl BioMart archive host so coordinates land on the build you actually - intend to intersect against The sc-linker enhancer-gene link files downloaded by - :func:`download_sclinker_enhancer_links` (Roadmap + ABC) are on - GRCh38, so pass ``genome_build="GRCh38"`` here (and to the LD/bim - panel used downstream) if you intend to keep everything on GRCh38 - instead of lifting the enhancer links over to GRCh37. + intend to intersect against. The sc-linker enhancer-gene link files + downloaded by :func:`download_sclinker_enhancer_links` (Roadmap + ABC) + are on **GRCh37** -- see :data:`SCLINKER_ENHANCER_LINKS_GENOME_BUILD` -- + so the default keeps everything on one build. Pass + ``genome_build="GRCh38"`` only if you lift the enhancer links over + first, and match the LD/bim panel used downstream. Requires ``pybiomart`` (``pip install pybiomart``). """ diff --git a/tests/test_sclinker_build_consistency.py b/tests/test_sclinker_build_consistency.py index 9d4c1f0..eb3fd0e 100644 --- a/tests/test_sclinker_build_consistency.py +++ b/tests/test_sclinker_build_consistency.py @@ -79,8 +79,9 @@ def test_genescores_to_100kb_bedgraph_propagates_genome_build(): assert bg.attrs.get("genome_build") == "GRCh37" -def test_sclinker_enhancer_links_build_is_grch38(): - assert SCLINKER_ENHANCER_LINKS_GENOME_BUILD == "GRCh38" +def test_sclinker_enhancer_links_build_is_grch37(): + """The Roadmap + ABC enhancer-gene links are distributed on GRCh37.""" + assert SCLINKER_ENHANCER_LINKS_GENOME_BUILD == "GRCh37" def test_bedgraph_to_snp_annotation_end_to_end_build_mismatch(tmp_path):