Source code for dae.gene_sets.gene_sets_db

"""Classes for handling of gene sets."""

import abc
import logging
import os
from functools import cached_property
from typing import Any, Literal

from pydantic import BaseModel, Field

from dae.gene_sets.gene_term import (
    read_ewa_set_file,
    read_gmt_file,
    read_mapping_file,
)
from dae.genomic_resources.fsspec_protocol import build_local_resource
from dae.genomic_resources.repository import (
    GenomicResource,
    GenomicResourceRepo,
)
from dae.genomic_resources.repository_factory import (
    build_genomic_resource_repository,
)
from dae.task_graph.graph import Task, TaskGraph

logger = logging.getLogger(__name__)


[docs] class MetaSchema(BaseModel): description: str | None = None labels: dict[str, Any] | None = None
[docs] class BaseResourceSchema(BaseModel): type: str | None = None meta: MetaSchema | None = None
[docs] class ViewRangeSchema(BaseModel): min: float | None = None max: float | None = None
# pylint: disable=missing-class-docstring
[docs] class HistogramSchema(BaseModel): type: str | None = None plot_function: str | None = None number_of_bins: int | None = None view_range: ViewRangeSchema | None = None x_log_scale: bool | None = None y_log_scale: bool | None = None x_min_log: float | None = None value_order: list[str | int] | None = None displayed_values_count: int | None = None displayed_values_percent: float | None = None reason: str | None = None
# pylint: disable=missing-class-docstring
[docs] class GeneSetResourceSchema(BaseResourceSchema): resource_id: str = Field(alias="id") filename: str | None = None directory: str | None = None resource_format: str | None = Field(alias="format") web_label: str | None = None web_format_str: str | None = None histograms: dict[ Literal["genes_per_gene_set", "gene_sets_per_gene"], HistogramSchema, ] | None = None
[docs] class GeneSet: """Class representing a set of genes.""" # pylint: disable=too-few-public-methods name: str desc: str count: int syms: list[str] def __init__(self, name: str, desc: str, syms: list[str]) -> None: self.name = name self.desc = desc self.count = len(syms) self.syms = syms def __getitem__(self, name: str) -> Any: # This is done so that GeneSet instances and # denovo gene set dictionaries can be accessed in a uniform way if name == "name": return self.name if name == "desc": return self.desc if name == "count": return self.count if name == "syms": return self.syms raise KeyError
[docs] class BaseGeneSetCollection(abc.ABC): """Base class for gene set collections."""
[docs] @abc.abstractmethod def get_gene_set(self, gene_set_id: str) -> GeneSet | None: """Return the gene set if found; returns None if not found.""" raise NotImplementedError
[docs] @abc.abstractmethod def get_all_gene_sets(self) -> list[GeneSet]: """Return list of all gene sets in the collection.""" raise NotImplementedError
[docs] class GeneSetCollection( BaseGeneSetCollection, ): """Class representing a collection of gene sets in a resource.""" def __init__(self, resource: GenomicResource) -> None: super().__init__() config = resource.get_config() if config is None: raise ValueError( f"genomic resource {resource.resource_id} not configured") self.resource = resource self.config = GeneSetResourceSchema.model_validate(config) self.collection_id = self.config.resource_id assert self.collection_id != "denovo" if resource.get_type() not in {"gene_set_collection", "gene_set"}: raise ValueError("Invalid resource type for gene set collection") if resource.get_type() == "gene_set": logger.warning( "'gene_set' resource type is deprecated; " "use 'gene_set_collection' instead") self.web_label = self.config.web_label self.web_format_str = self.config.web_format_str logger.debug("loading %s: %s", self.collection_id, config) self.gene_sets: dict[str, GeneSet] = self.load_gene_sets() assert self.collection_id, self.gene_sets @property def files(self) -> set[str]: """Return a list of resource files the implementation utilises.""" res = set() collection_format = self.config.resource_format if collection_format == "map": filename = self.config.filename assert filename is not None res.add(filename) names_filename = filename[:-4] + "names.txt" if self.resource.file_exists(names_filename): res.add(names_filename) elif collection_format == "gmt": filename = self.config.filename assert filename is not None res.add(filename) elif collection_format == "directory": directory = self.config.directory assert directory is not None if directory == ".": directory = "" for filepath, _ in self.resource.get_manifest().get_files(): if filepath.startswith(directory) and \ filepath.endswith(".txt"): res.add(filepath) else: raise ValueError("Invalid collection format type") return res
[docs] def load_gene_sets(self) -> dict[str, GeneSet]: """Build a gene set collection from a given GenomicResource.""" assert self.resource is not None gene_sets = {} collection_format = self.config.resource_format logger.debug("loading %s", self.collection_id) if collection_format == "map": filename = self.config.filename assert filename is not None names_filename = filename[:-4] + "names.txt" names_file = None if self.resource.file_exists(names_filename): names_file = self.resource.open_raw_file(names_filename) gene_terms = read_mapping_file( self.resource.open_raw_file(filename), names_file, ) elif collection_format == "gmt": filename = self.config.filename assert filename is not None gene_terms = read_gmt_file(self.resource.open_raw_file(filename)) elif collection_format == "directory": directory = self.config.directory assert directory is not None filepaths = [] if directory == ".": directory = "" # Easier check with startswith for filepath, _ in self.resource.get_manifest().get_files(): if filepath.startswith(directory) and \ filepath.endswith(".txt"): filepaths.append(filepath) files = [self.resource.open_raw_file(f) for f in filepaths] gene_terms = read_ewa_set_file(files) else: raise ValueError("Invalid collection format type") for key, value in gene_terms.t_desc.items(): syms = list(gene_terms.t2g[key].keys()) gene_set = GeneSet(key, value, syms) gene_sets[gene_set.name] = gene_set return gene_sets
[docs] def get_gene_set(self, gene_set_id: str) -> GeneSet | None: """Return the gene set if found; returns None if not found.""" gene_set = self.gene_sets.get(gene_set_id) if gene_set is None: logger.warning( "%s not found in %s", gene_set_id, self.gene_sets.keys(), ) return gene_set
[docs] def get_all_gene_sets(self) -> list[GeneSet]: return list(self.gene_sets.values())
[docs] def calc_info_hash(self) -> bytes: return b"placeholder"
[docs] def calc_statistics_hash(self) -> bytes: return b"placeholder"
[docs] def add_statistics_build_tasks( self, task_graph: TaskGraph, **kwargs: Any, # noqa: ARG002 ) -> list[Task]: return []
[docs] def get_genes_per_gene_set_hist_filename(self) -> str: return "statistics/genes_per_gene_set_histogram.png"
[docs] def get_gene_sets_per_gene_hist_filename(self) -> str: return "statistics/gene_sets_per_gene_histogram.png"
[docs] class GeneSetsDb: """Class that represents a dictionary of gene set collections.""" def __init__(self, gene_set_collections: list[GeneSetCollection]) -> None: self.gene_set_collections: dict[str, GeneSetCollection] = { gsc.collection_id: gsc for gsc in gene_set_collections } @cached_property def collections_descriptions(self) -> list[dict[str, Any]]: """Collect gene set descriptions. Iterates and creates a list of descriptions for each gene set collection """ gene_sets_collections_desc = [] for gsc in self.gene_set_collections.values(): label = gsc.web_label format_str = gsc.web_format_str gsc_id = gsc.collection_id if not label or not format_str: continue gene_sets_collections_desc.append( { "desc": label, "name": gsc_id, "format": format_str.split("|"), }, ) return gene_sets_collections_desc
[docs] def has_gene_set_collection(self, gsc_id: str) -> bool: """Check the database if contains the specified gene set collection.""" return gsc_id in self.gene_set_collections
[docs] def get_gene_set_collection_ids(self) -> set[str]: """Return all gene set collection ids. Including the ids of collections which have not been loaded. """ return set(self.gene_set_collections.keys())
[docs] def get_gene_set_ids(self, collection_id: str) -> set[str]: """Return the IDs of all the gene sets in specified collection.""" gsc = self.gene_set_collections[collection_id] return set(gsc.gene_sets.keys())
[docs] def get_all_gene_sets(self, collection_id: str) -> list[GeneSet]: """Return all the gene sets in the specified collection.""" gsc = self.gene_set_collections[collection_id] logger.debug( "gene sets from %s: %s", collection_id, len(gsc.gene_sets.keys())) return list(gsc.gene_sets.values())
[docs] def get_gene_set( self, collection_id: str, gene_set_id: str, ) -> GeneSet | None: """Find and return a gene set in a gene set collection.""" gsc = self.gene_set_collections[collection_id] return gsc.get_gene_set(gene_set_id)
def __len__(self) -> int: return len(self.gene_set_collections)
[docs] def build_gene_set_collection_from_file( filename: str, collection_id: str | None = None, collection_format: str | None = None, web_label: str | None = None, web_format_str: str | None = None, ) -> GeneSetCollection: """Return a Gene Set Collection by adapting a file to a local resource.""" dirname = os.path.dirname(filename) basename = os.path.basename(filename) if collection_format is None: is_dir = os.path.isdir(filename) if is_dir: collection_format = "directory" else: extension = os.path.splitext(filename)[1] if extension == ".txt": collection_format = "map" elif extension == ".gmt": collection_format = "gmt" elif extension == ".sql": collection_format = "sqlite" else: raise ValueError("Cannot find collection format automatically") if collection_id is None: collection_id = basename config = { "type": "gene_set", "id": collection_id, "format": collection_format, "web_label": web_label, "web_format_str": web_format_str, } if collection_format == "directory": config["directory"] = basename elif collection_format == "sqlite": config["dbfile"] = basename else: config["filename"] = basename resource = build_local_resource(dirname, config) return build_gene_set_collection_from_resource(resource)
[docs] def build_gene_set_collection_from_resource( resource: GenomicResource, ) -> GeneSetCollection: """Return a Gene Set Collection built from a resource.""" if resource is None: raise ValueError(f"missing resource {resource}") return GeneSetCollection(resource)
[docs] def build_gene_set_collection_from_resource_id( resource_id: str, grr: GenomicResourceRepo | None = None, ) -> GeneSetCollection: if grr is None: grr = build_genomic_resource_repository() return build_gene_set_collection_from_resource( grr.get_resource(resource_id))