diff --git a/pyproject.toml b/pyproject.toml index e84dec9..24b188f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,12 +30,16 @@ classifiers = [ "Topic :: Scientific/Engineering :: Bio-Informatics", ] dependencies = [ + "anndata>=0.11.4", + "boltons>=25.0.0", # ctxcore dependency + "ctxcore", "ncls>=0.0.70", "polars>=1.31.0", "pyarrow>=20.0.0", "pybiomart>=0.2.0", "requests>=2.32.4", "scipy>=1.15.3", + "setuptools<81", # ctxcore dependency "tmtoolkit>=0.12.0", ] @@ -126,3 +130,6 @@ ban-relative-imports = "all" [tool.mypy] mypy_path = "$MYPY_CONFIG_FILE_DIR/stubs" plugins = ["numpy.typing.mypy_plugin"] + +[tool.uv.sources] +ctxcore = { git = "https://github.com/aertslab/ctxcore.git" } diff --git a/src/pycisTopic/cli/pycistopic.py b/src/pycisTopic/cli/pycistopic.py index 37f4c9b..fa5068b 100644 --- a/src/pycisTopic/cli/pycistopic.py +++ b/src/pycisTopic/cli/pycistopic.py @@ -4,6 +4,9 @@ from pycisTopic.cli.subcommand.count_matrix import add_parser_count_matrix from pycisTopic.cli.subcommand.qc import add_parser_qc +from pycisTopic.cli.subcommand.signature_enrichment import ( + add_parser_signature_enrichment, +) from pycisTopic.cli.subcommand.topic_modeling import add_parser_topic_modeling from pycisTopic.cli.subcommand.tss import add_parser_tss @@ -23,6 +26,7 @@ def main(): add_parser_topic_modeling(subparsers) add_parser_tss(subparsers) add_parser_count_matrix(subparsers) + add_parser_signature_enrichment(subparsers) args = parser.parse_args() args.func(args) diff --git a/src/pycisTopic/cli/subcommand/signature_enrichment.py b/src/pycisTopic/cli/subcommand/signature_enrichment.py new file mode 100644 index 0000000..4f028e4 --- /dev/null +++ b/src/pycisTopic/cli/subcommand/signature_enrichment.py @@ -0,0 +1,220 @@ +from argparse import _SubParsersAction + +import polars as pl + +OBSM_SIGN_ENRICHMENT_KEY="AUCell" +UNS_SIGN_ENRICHMENT_KEY="AUCell signature names" + + +def _region_names_to_polars_df( + region_names: list[str], + chrom_split: str = ":", + start_end_split: str = "-", +) -> pl.DataFrame: + chromosomes: list[str] = [] + starts: list[int] = [] + ends: list[int] = [] + for region in region_names: + chrom, start, end = region.replace( + chrom_split, start_end_split + ).split(start_end_split) + chromosomes.append(chrom) + starts.append(int(start)) + ends.append(int(end)) + + return pl.from_dict( + { + "Chromosome": chromosomes, + "Start": starts, + "End": ends + } + ) + +def run_signature_enrichment(args): + import logging + import sys + + import anndata as ad # type: ignore + import numpy as np + + from pycisTopic.fragments import read_bed_to_polars_df + from pycisTopic.signature_enrichment import signature_enrichment + + if len(args.signature_names) != len(args.signatures): + raise ValueError( + f"Length of signature_names ({len(args.signature_names)}) does not " + f"match the length of the signatures ({len(args.signatures)})!" + ) + + if args.verbose: + level = logging.INFO + log_format = "%(asctime)s %(name)-12s %(levelname)-8s %(message)s" + handlers = [logging.StreamHandler(stream=sys.stdout)] + logging.basicConfig(level=level, format=log_format, handlers=handlers) + log = logging.getLogger("cisTopic") + else: + log = None + + if log is not None: + log.info("Reading data ...") + + ad_region_topic = ad.read_h5ad(args.region_topic) + ad_cell_topic = ad.read_h5ad(args.cell_topic) + + assert isinstance(ad_region_topic.X, np.ndarray), \ + "`ad_region_topic.X` should be numpy array" + + assert isinstance(ad_cell_topic.X, np.ndarray), \ + "`ad_cell_topic.X` should be numpy array" + + region_topic_granges = _region_names_to_polars_df( + ad_region_topic.obs_names # type: ignore + ) + + signatures: dict[str, pl.DataFrame] = {} + for signature_name, signature_bed in zip( + args.signature_names, args.signatures, + strict=True + ): + signatures[signature_name] = read_bed_to_polars_df(signature_bed) + + auc_values = signature_enrichment( + region_topic=ad_region_topic.X, + cell_topic=ad_cell_topic.X.T, + region_topic_granges=region_topic_granges, + signatures=signatures, + chunk_size=args.chunk_size, + normalize=args.normalize, + min_frac_consensus=args.min_frac_consensus, + min_frac_signature=args.min_frac_signature, + seed=args.seed, + auc_threshold=args.auc_threshold, + log=log + ) + + if OBSM_SIGN_ENRICHMENT_KEY not in ad_cell_topic.obsm: + if log is not None: + log.info( + f"Storing AUCell matrix in .obsm[{OBSM_SIGN_ENRICHMENT_KEY}] " + f"and signature names in .uns[{UNS_SIGN_ENRICHMENT_KEY}]." + ) + ad_cell_topic.obsm[OBSM_SIGN_ENRICHMENT_KEY] = auc_values + ad_cell_topic.uns[UNS_SIGN_ENRICHMENT_KEY] = list(signatures.keys()) + else: + if log is not None: + log.info( + f"Concatenating new AUCell matrix with .obsm[{OBSM_SIGN_ENRICHMENT_KEY}] " + f"and extending signature names in .uns[{UNS_SIGN_ENRICHMENT_KEY}]." + ) + ad_cell_topic.obsm[OBSM_SIGN_ENRICHMENT_KEY] = np.concatenate( + (ad_cell_topic.obsm[OBSM_SIGN_ENRICHMENT_KEY], auc_values), + axis=1 + ) + assert isinstance(ad_cell_topic.uns[UNS_SIGN_ENRICHMENT_KEY], list), \ + "cell topic anndata contains AUCell values but no signature names in " + \ + f".uns[{UNS_SIGN_ENRICHMENT_KEY}]." + ad_cell_topic.uns[UNS_SIGN_ENRICHMENT_KEY].extend(list(signatures.keys())) + + if log is not None: + log.info(f"Saving cell topic anndata to: {args.out_file}") + ad_cell_topic.write(args.out_file) + +def add_parser_signature_enrichment( + subparser: _SubParsersAction +): + parser_signature_enrichment = subparser.add_parser( + "signature_enrichment", + help="Run signature enrichment using AUCell." + ) + parser_signature_enrichment.set_defaults(func=run_signature_enrichment) + + parser_signature_enrichment.add_argument( + "-r", + "--region-topic", + dest="region_topic", + type=str, + required=True, + help="Path to region-topic AnnData (.h5ad) file." + ) + parser_signature_enrichment.add_argument( + "-c", + "--cell-topic", + dest="cell_topic", + type=str, + required=True, + help="Path to cell-topic AnnData (.h5ad) file." + ) + parser_signature_enrichment.add_argument( + "-n", + "--signature-names", + dest="signature_names", + nargs="+", + type=str, + required=True, + help="List of signature names (space-separated, same order as --signatures)." + ) + parser_signature_enrichment.add_argument( + "-s", + "--signatures", + dest="signatures", + nargs="+", + type=str, + required=True, + help="List of BED files for signatures (space-separated, same order as --signature-names)." + ) + parser_signature_enrichment.add_argument( + "-o", + "--out-file", + dest="out_file", + type=str, + required=True, + help="Output file for updated cell-topic AnnData (.h5ad)." + ) + parser_signature_enrichment.add_argument( + "--chunk-size", + dest="chunk_size", + type=int, + default=1000, + help="Number of cells to process at once. Default: 1000." + ) + parser_signature_enrichment.add_argument( + "--normalize", + dest="normalize", + action="store_true", + help="Normalize AUC values to a maximum of 1.0 per regulon." + ) + parser_signature_enrichment.add_argument( + "--min-frac-consensus", + dest="min_frac_consensus", + type=float, + default=0.5, + help="Minimal fractional overlap of signature and consensus peak relative to consensus peak. Default: 0.5." + ) + parser_signature_enrichment.add_argument( + "--min-frac-signature", + dest="min_frac_signature", + type=float, + default=0.5, + help="Minimal fractional overlap of signature and consensus peak relative to signature. Default: 0.5." + ) + parser_signature_enrichment.add_argument( + "--seed", + dest="seed", + type=int, + default=12345, + help="Random seed for ranking. Default: 12345." + ) + parser_signature_enrichment.add_argument( + "--auc-threshold", + dest="auc_threshold", + type=float, + default=0.05, + help="Fraction of ranked genome for AUC calculation. Default: 0.05." + ) + parser_signature_enrichment.add_argument( + "-v", + "--verbose", + dest="verbose", + action="store_true", + help="Enable verbose logging." + ) diff --git a/src/pycisTopic/imputed_accessibility.py b/src/pycisTopic/imputed_accessibility.py new file mode 100644 index 0000000..db0bf2c --- /dev/null +++ b/src/pycisTopic/imputed_accessibility.py @@ -0,0 +1,242 @@ +import logging +from typing import Iterator, Literal + +import numpy as np +import numpy.typing as npt +import polars as pl + + +def _alocate_chunk_array( + n_cells: int, + n_regions: int, + chunk_size: int, + chunk_along: Literal["cell", "region"], + log: logging.Logger | None +) -> npt.NDArray[np.float32]: + if chunk_along == "region": + if log is not None: + log.info( + f"Allocate {(chunk_size * n_cells * 4 / 1024**3):.3f} GiB of RAM for " + f"calculating (partial) imputed accessibility per cell for ({n_cells}) cells " + f"for chunk of {chunk_size} regions." + ) + return np.empty( + (chunk_size, n_cells), + dtype=np.float32 + ) + else: + if log is not None: + log.info( + f"Allocate {(chunk_size * n_regions * 4 / 1024**3):.3f} GiB of RAM for " + f"calculating (partial) imputed accessibility per region for ({n_regions}) regions " + f"for chunk of {chunk_size} cells." + ) + return np.empty( + (n_regions, chunk_size), + dtype=np.float32 + ) + +def impute_accessibility_chunked( + region_topic: npt.NDArray[np.float32], + cell_topic: npt.NDArray[np.float32], + chunk_size: int, + chunk_along: Literal["cell", "region"], + log: logging.Logger | None = None, + return_start_end: bool = False, +) -> Iterator[npt.NDArray[np.float32]] | Iterator[tuple[tuple[int, int], npt.NDArray[np.float32]]]: + """ + Impute accessibility in chunks. + + Parameters + ---------- + region_topic + Region topic matrix (regions x topics). + cell_topic + Cell topic matrix (topic x cells). + chunk_size + The size of the chunks. + chunk_along + Whether to chunk along cells or regions. + log + Optional logging.Logger + return_start_end + Whether to return the start and end index of the chunk. + + Yields + ------ + Numpy array with imputed accessibility of chunk. + + .. warning:: + The yielded array is a reused in-place buffer. Each call to ``next()`` + overwrites its contents. Do not hold references to previously yielded + arrays across iterations — copy the array (e.g. ``chunk.copy()``) if + you need to retain the data. + + """ + if region_topic.shape[1] != cell_topic.shape[0]: + raise ValueError( + "region- and cell-topic dimensions do not match" + f" region_topic: {region_topic.shape[1]}" + f" cell_topic: {cell_topic.shape[0]}." + "Maybe you have to transpose region_topic or cell_topic" + ) + + if chunk_along not in ["cell", "region"]: + raise ValueError(f"chunk along should be either 'cell' or 'region', not {chunk_along}!") + + if log is not None: + log.info(f"Calculating imputed accessibility in chunks across {chunk_along} ...") + + n_regions = region_topic.shape[0] + n_cells = cell_topic.shape[1] + + imputed_acc_chunk: npt.NDArray[np.float32] = _alocate_chunk_array( + n_cells=n_cells, + n_regions=n_regions, + chunk_size=chunk_size, + chunk_along=chunk_along, + log=log + ) + + n_total: int = n_regions if chunk_along == "region" else n_cells + + for chunk_start in range(0, n_total, chunk_size): + chunk_end = chunk_start + chunk_size + + # get current chunk of regions or cells + if chunk_along == "region": + _cell_topic = cell_topic + _region_topic = region_topic[ + chunk_start: chunk_end + ] + current_chunk_size = _region_topic.shape[0] + else: + _cell_topic = cell_topic[:, chunk_start: chunk_end] + _region_topic = region_topic + current_chunk_size = _cell_topic.shape[1] + + if current_chunk_size < chunk_size: + del imputed_acc_chunk + imputed_acc_chunk = _alocate_chunk_array( + n_cells=n_cells, + n_regions=n_regions, + chunk_size=current_chunk_size, + chunk_along=chunk_along, + log=log + ) + chunk_end = chunk_start + current_chunk_size + + if log is not None: + log.info( + "Calculate partial imputed accessibility " + f"{chunk_start}-{chunk_end} (out of {n_total})." + ) + + np.matmul(_region_topic, _cell_topic, out=imputed_acc_chunk) + if return_start_end: + yield (chunk_start, chunk_end), imputed_acc_chunk + else: + yield imputed_acc_chunk + +def rank_imputed_accessibility( + imputed_accessibility: npt.NDArray[np.float32], + seed: int = 123, + method: Literal["numpy", "polars"] = "polars", +) -> npt.NDArray[np.int32]: + """ + Generate rankings per cell based on the imputed accessibility scores per region. + + Parameters + ---------- + imputed_accessibility + Numpy array of imputed accessibility (regions x cells) + seed + Random seed to ensure reproducibility of the rankings when there are ties + method + Method to use for the ranking implementation. + Options are "numpy" (same rankings as with older versions of pycisTopic) + or "polars" (fastest). Default: "polars". + + Return + ------ + Numpy array + Containing ranking values rather than scores (regions x cells). + + """ + if method != "numpy" and method != "polars": + raise ValueError( + f'Invalid method ("{method}") for ranking implementation. Use "numpy" or "polars".' + ) + + # Initialize random number generator, for handling ties. + rng = np.random.default_rng(seed=seed) + + # Function to make rankings per array. + def rank_scores_and_assign_random_ranking_in_range_for_ties_with_numpy( + scores_with_ties_for_motif_or_track_numpy: npt.NDArray[np.float32], + ) -> npt.NDArray[np.int32]: + # + # Create random permutation so tied scores will have a different ranking each time. + random_permutations_to_break_ties_numpy = rng.permutation( + scores_with_ties_for_motif_or_track_numpy.shape[0] + ) + ranking_with_broken_ties_for_motif_or_track_numpy = np.empty( + scores_with_ties_for_motif_or_track_numpy.shape[0], + dtype=np.int32, + ) + ranking_with_broken_ties_for_motif_or_track_numpy[ + random_permutations_to_break_ties_numpy[ + (-scores_with_ties_for_motif_or_track_numpy)[ + random_permutations_to_break_ties_numpy + ].argsort() + ] + ] = np.arange( + scores_with_ties_for_motif_or_track_numpy.shape[0], + dtype=np.int32, + ) + + return ranking_with_broken_ties_for_motif_or_track_numpy + + def rank_scores_and_assign_random_ranking_in_range_for_ties_with_polars( + scores_with_ties_for_motif_or_track_numpy: npt.NDArray[np.float32], + seed: int, + ) -> npt.NDArray[np.int32]: + # Rank scores and assign a random ranking in range for regions/genes with + # the same score. + # - Convert numpy array to Polars Series . + # - Replace NaN values with the minimum value of the dtype, so NaNs are + # ranked last. + # - Use the `rank` method from Polars to assign ranks, using the "random" + # method to break ties so that regions/genes with the same score get a + # random ranking in the range of their scores instead of depending on + # the order in which they appear in the input array. + # - Subtract 1 from the ranks to make them zero-based. + return ( + pl.Series(scores_with_ties_for_motif_or_track_numpy) + .fill_nan(float(np.finfo(scores_with_ties_for_motif_or_track_numpy.dtype).min)) + .rank(method="random", descending=True, seed=seed) + - 1 + ).to_numpy() + + n_regions, n_cells = imputed_accessibility.shape + # Create zeroed imputed object rankings database. + rankings = np.zeros((n_regions, n_cells), dtype=np.int32) + + # Rank all scores per motif/track and assign a random ranking in range for regions/genes with the same score. + if method == "numpy": + for cell_idx in range(n_cells): + rankings[:, cell_idx] = ( + rank_scores_and_assign_random_ranking_in_range_for_ties_with_numpy( + imputed_accessibility[:, cell_idx] + ) + ) + else: + for cell_idx in range(n_cells): + rankings[:, cell_idx] = ( + rank_scores_and_assign_random_ranking_in_range_for_ties_with_polars( + imputed_accessibility[:, cell_idx], + seed=rng.integers(2**32 - 1), + ) + ) + + return rankings diff --git a/src/pycisTopic/signature_enrichment.py b/src/pycisTopic/signature_enrichment.py index 28a7314..de4a536 100644 --- a/src/pycisTopic/signature_enrichment.py +++ b/src/pycisTopic/signature_enrichment.py @@ -1,136 +1,179 @@ from __future__ import annotations +from multiprocessing import cpu_count from typing import TYPE_CHECKING -import pyranges as pr -from ctxcore.genesig import GeneSignature -from pyscenic.aucell import aucell4r +import numpy as np +import numpy.typing as npt +import pandas as pd # type: ignore +import polars as pl +from ctxcore.aucell import aucell4r # type: ignore +from ctxcore.genesig import GeneSignature # type: ignore -if TYPE_CHECKING: - from pycisTopic.diff_features import CistopicImputedFeatures - -# FIXME -from .diff_features import * -from .utils import * +from pycisTopic.genomic_ranges import intersection +from pycisTopic.imputed_accessibility import ( + impute_accessibility_chunked, + rank_imputed_accessibility, +) +if TYPE_CHECKING: + import logging + + +def _polars_granges_to_region_names( + granges: pl.DataFrame, +) -> list[str]: + return granges.with_columns( + region_names = ( + pl.concat_str( + [ + pl.col("Chromosome"), + pl.concat_str( + [ + pl.col("Start"), + pl.col("End") + ], + separator="-" + ) + ], + separator=":" + ) + ) + )["region_names"].to_list() + +def _region_names_to_signature( + region_names: list[str], + name: str, +) -> GeneSignature: + weights = np.ones(len(region_names)) + return GeneSignature( + name=name, gene2weight=dict(zip(region_names, weights)) + ) def signature_enrichment( - rankings: CistopicImputedFeatures, - signatures: dict[str, pr.PyRanges] | dict[str, list], - enrichment_type: str = "region", + region_topic: npt.NDArray[np.float32], + cell_topic: npt.NDArray[np.float32], + region_topic_granges: pl.DataFrame, + signatures: dict[str, pl.DataFrame], + chunk_size: int, + normalize: bool, + min_frac_consensus: float, + min_frac_signature: float, + seed: int, auc_threshold: float = 0.05, - normalize: bool = False, - n_cpu: int = 1, -): + log: logging.Logger | None = None, +) -> npt.NDArray[np.float32]: """ - Get enrichment of a region signature in cells or topics using AUCell (Van de Sande et al., 2020). + Calculate enrichment of region signatures in cells using AUCell (Van de Sande et al., 2020). Parameters ---------- - rankings: CistopicImputedFeatures - A CistopicImputedFeatures object with ranking values - signatures: Dictionary of pr.PyRanges (for regions) or list (for genes) - A dictionary containing region signatures as pr.PyRanges or gene names as list - enrichment_type: str - Whether features are genes or regions - auc_threshold: float - The fraction of the ranked genome to take into account for the calculation of the Area Under the recovery Curve. Default: 0.05 + region_topic + Region topic matrix (regions x topics). + cell_topic + Cell topic matrix (topic x cells). + region_topic_granges + Polars Dataframe with genomic ranges corresponding to region topic. + signatures + Dictionary of genomic ranges signatures (polars DataFrames). + chunk_size + The number of cells to process at once. normalize: bool - Normalize the AUC values to a maximum of 1.0 per regulon. Default: False - num_workers: int - The number of cores to use. Default: 1 - - Return - ------ - A pd.DataFrame containing signatures as columns, cells/topics as rows and AUC scores as values + Normalize the AUC values to a maximum of 1.0 per regulon. + min_frac_consensus + Minimal fractional overlap of signature and consensus peak + relative to consensus peak. + min_frac_signature + Minimal fractional overlap of signature and consensus peak + relative to signature. + seed + Seed used to randomly resolve tied values in imputed accessibility + for generatin the ranking. + auc_threshold: float + The fraction of the ranked genome to take into account for the calculation + of the Area Under the recovery Curve. Default: 0.05. + log + Optional logger. - References - ---------- - Van de Sande, B., Flerin, C., Davie, K., De Waegeneer, M., Hulselmans, G., Aibar, S., ... & Aerts, S. (2020). A scalable SCENIC workflow for single-cell gene - regulatory network analysis. Nature Protocols, 15(7), 2247-2276. + Returns + ------- + A Numpy array with auc values across cells (cell x signature) """ - # Compute rankings if needed and format input - rankings = pd.DataFrame( - rankings.mtx.transpose(), - columns=rankings.feature_names, - index=rankings.cell_names, - ) - # Take regions in input - if enrichment_type == "region": - regions_in_input = pr.PyRanges( - region_names_to_coordinates(rankings.columns.tolist()) + if region_topic_granges.shape[0] != region_topic.shape[0]: + raise ValueError( + f"Length of the region names ({region_topic_granges.shape[0]}) " + f"does not match the shape of region_topic {region_topic.shape}" + ) + region_names = _polars_granges_to_region_names(region_topic_granges) + + # Put signatures in coordinate frame of region topic by performing + # intersect and retaining regions of the region topic granges that pass + # the overlap thresholds (min_frac_consensus and min_frac_signatures). + gr_signatures_consensus: dict[str, pl.DataFrame] = {} + for signature, sign_granges in signatures.items(): + gr_signatures_consensus[signature] = intersection( + regions1_df_pl=region_topic_granges, + regions2_df_pl=sign_granges, + regions1_coord=True, + add_overlap_size=True, + regions1_suffix="@1" + ).filter( + (pl.col("fr1_inter") >= min_frac_consensus) & + (pl.col("fr2_inter") >= min_frac_signature) + ).select( + pl.col("Chromosome@1"), + pl.col("Start@1"), + pl.col("End@1") + ).rename( + { + "Chromosome@1": "Chromosome", + "Start@1": "Start", + "End@1": "End" + } ) - # Get signatures - signatures = [ - region_set_to_signature(signatures[key], regions_in_input, key) - for key in signatures.keys() - ] - if enrichment_type == "gene": - # Get signatures - signatures = [ - gene_set_to_signature(signatures[key], key) for key in signatures.keys() - ] - # Run aucell - auc_sig = aucell4r( - df_rnk=rankings, - signatures=signatures, - auc_threshold=auc_threshold, - noweights=False, - normalize=normalize, - num_workers=n_cpu, - ) - auc_sig.columns.names = [None] - return auc_sig - - -def region_set_to_signature( - query_region_set: pr.PyRanges, target_region_set: pr.PyRanges, name: str -): - """ - A helper function to intersect query regions with the input data set regions. - - Parameters - ---------- - query_region_set: pr.PyRanges - Pyranges with regions to query - target_region_set: pr.PyRanges - Pyranges with target regions - name: str - Name for the signature - - Return - ------ - A GeneSignature object to use with AUCell - """ - query_in_target = query_region_set.join(target_region_set) - query_in_target = query_in_target.df[["Chromosome", "Start_b", "End_b"]] - query_in_target.columns = ["Chromosome", "Start", "End"] - query_in_target = coord_to_region_names(pr.PyRanges(query_in_target)) - weights = np.ones(len(query_in_target)) - signature = GeneSignature( - name=name, gene2weight=dict(zip(query_in_target, weights)) + # Convert granges signatures to gene signatures + gs_signatures_consensus: list[GeneSignature] = [ + _region_names_to_signature( + region_names=_polars_granges_to_region_names(granges), + name=signature + ) + for signature, granges in gr_signatures_consensus.items() + ] + + # initialize aucell values + n_cells = cell_topic.shape[1] + n_signatures = len(signatures) + aucell_values: npt.NDArray[np.float32] = np.empty( + (n_cells, n_signatures), dtype=np.float32 ) - return signature - - -def gene_set_to_signature(gene_set: list, name: str): - """ - A helper function to generat gene signatures. - - Parameters - ---------- - gene_set: pr.PyRanges - List of genes - name: str - Name for the signature - - Return - ------ - A GeneSignature object to use with AUCell - - """ - weights = np.ones(len(gene_set)) - signature = GeneSignature(name=name, gene2weight=dict(zip(gene_set, weights))) - return signature + for (cell_start, cell_end), imputed_acc_chunk in impute_accessibility_chunked( + region_topic=region_topic, + cell_topic=cell_topic, + chunk_size=chunk_size, + chunk_along="cell", + log=log, + return_start_end=True, + ): + if log is not None: + log.info("Generating ranking.") + ranking_chunk = rank_imputed_accessibility( + imputed_accessibility=imputed_acc_chunk, + seed=seed + ) + if log is not None: + log.info("Calculating AUCs.") + aucell_values[cell_start: cell_end] = aucell4r( + df_rnk=pd.DataFrame( + ranking_chunk.T, + columns=region_names + ), + signatures=gs_signatures_consensus, # type: ignore + auc_threshold=auc_threshold, + noweights=False, + normalize=False, + num_workers=min(chunk_size, cpu_count()), + ).to_numpy() + + return aucell_values / aucell_values.max(axis=0) if normalize else aucell_values