Genomic analysis

Driven adapters that implement genomics domain ports backed by third-party scientific libraries (BiocPy, pyhpo), plus concrete GenotypeProcessor implementations for cohort-level statistics injected into ZarrVariantStore.

Overview

Ports implemented

Adapter	Port	Backend
BiocpyRangeOperations	RangeOperations	BiocPy GenomicRanges
PyHpoOperations	HpoOperations	pyhpo

BiocpyRangeOperations

Bridges GenomicPosition domain objects and BiocPy's GenomicRanges. All coordinate conversion is handled internally: iris is 1-based inclusive, BiocPy GenomicRanges is 0-based half-open (iris start=1, end=10 ↔ BiocPy start=0, end=10).

PyHpoOperations

Implements HpoOperations over the pyhpo ontology. Term resolution prefers hpo_id when available, falling back to name otherwise. pyhpo initializes a module-level singleton on the first Ontology() call — subsequent PyHpoOperations instances reuse the same in-memory ontology rather than re-reading it from disk. data_directory may be passed to point at a custom HPO data dump; when None (default), pyhpo's bundled data is used.

GenotypeProcessor implementations

Concrete implementations of the GenotypeProcessor port, passed to ZarrVariantStore(processor=...) — see Stores & Pipeline.

Processor	Adds to VariantAnnotation.metadata
AncestryStratifiedAFProcessor	Per-decile allele frequency stratified by an ancestry covariate, plus a linear regression (slope, SE, p-value, R²) of AF against ancestry proportion
PerIndividualAltCountProcessor	Median and IQR of alt-allele counts and carrier status per individual across a variant block

BiocPy range operations

iris.adapters.driven.genomic_analysis.biocpy

BiocPy adapter for the RangeOperations port.

Bridges between iris GenomicPosition domain objects and the BiocPy GenomicRanges library. All coordinate conversion (1-based inclusive ↔ 0-based half-open) is handled here, keeping the domain model clean.

BiocPy GenomicRanges uses 0-based half-open intervals internally

• iris start=1, end=10 → BiocPy start=0, end=10
• BiocPy start=0, end=10 → iris start=1, end=10

BiocpyRangeOperations

BiocpyRangeOperations(default_assembly=None)

Bases: RangeOperations

RangeOperations adapter backed by BiocPy GenomicRanges.

Converts between iris GenomicPosition objects and BiocPy GenomicRanges, delegating all computation to the BiocPy library.

PARAMETER	DESCRIPTION
default_assembly	Assembly name to attach to returned GenomicPosition objects. If None, assembly is not set. TYPE: str | None DEFAULT: None

Example

from iris.domain.genomics.services import GenomicRangeService

ops = BiocpyRangeOperations()
service = GenomicRangeService(ops=ops)
merged = service.union(query_ranges, subject_ranges)

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def __init__(self, default_assembly: str | None = None) -> None:
    self._assembly = default_assembly

union

union(a, b)

Return the union of two sets of genomic ranges.

PARAMETER	DESCRIPTION
a	First set of genomic ranges. TYPE: Sequence[GenomicPosition]
b	Second set of genomic ranges. TYPE: Sequence[GenomicPosition]

RETURNS	DESCRIPTION
list[GenomicPosition]	Merged, non-overlapping ranges covering all positions in a or b.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def union(
    self, a: Sequence[GenomicPosition], b: Sequence[GenomicPosition]
) -> list[GenomicPosition]:
    """Return the union of two sets of genomic ranges.

    Args:
        a: First set of genomic ranges.
        b: Second set of genomic ranges.

    Returns:
        Merged, non-overlapping ranges covering all positions in a or b.
    """
    return self._pos(self._gr(a).union(self._gr(b)))

intersect

intersect(a, b)

Return ranges present in both a and b.

PARAMETER	DESCRIPTION
a	First set of genomic ranges. TYPE: Sequence[GenomicPosition]
b	Second set of genomic ranges. TYPE: Sequence[GenomicPosition]

RETURNS	DESCRIPTION
list[GenomicPosition]	Ranges representing the intersection of a and b.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def intersect(
    self, a: Sequence[GenomicPosition], b: Sequence[GenomicPosition]
) -> list[GenomicPosition]:
    """Return ranges present in both a and b.

    Args:
        a: First set of genomic ranges.
        b: Second set of genomic ranges.

    Returns:
        Ranges representing the intersection of a and b.
    """
    return self._pos(self._gr(a).intersect(self._gr(b)))

setdiff

setdiff(a, b)

Return ranges in a not covered by b.

PARAMETER	DESCRIPTION
a	Query ranges. TYPE: Sequence[GenomicPosition]
b	Ranges to subtract. TYPE: Sequence[GenomicPosition]

RETURNS	DESCRIPTION
list[GenomicPosition]	Ranges in a not overlapping any range in b.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def setdiff(
    self, a: Sequence[GenomicPosition], b: Sequence[GenomicPosition]
) -> list[GenomicPosition]:
    """Return ranges in a not covered by b.

    Args:
        a: Query ranges.
        b: Ranges to subtract.

    Returns:
        Ranges in a not overlapping any range in b.
    """
    return self._pos(self._gr(a).setdiff(self._gr(b)))

coverage

coverage(ranges)

Compute per-base coverage across all ranges.

PARAMETER	DESCRIPTION
ranges	Genomic ranges to compute coverage for. TYPE: Sequence[GenomicPosition]

RETURNS	DESCRIPTION
dict[str, list[int]]	Dict mapping chromosome to list of per-base coverage values.
dict[str, list[int]]	Index 0 corresponds to position 1 on that chromosome.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def coverage(self, ranges: Sequence[GenomicPosition]) -> dict[str, list[int]]:
    """Compute per-base coverage across all ranges.

    Args:
        ranges: Genomic ranges to compute coverage for.

    Returns:
        Dict mapping chromosome to list of per-base coverage values.
        Index 0 corresponds to position 1 on that chromosome.
    """
    gr = self._gr(ranges)
    cov = gr.coverage()
    # BiocPy returns a RleList-like object; convert to plain dict of lists
    return {chrom: list(cov[chrom]) for chrom in cov.get_names()}

gaps

gaps(ranges, *, genome=None)

Return gaps between ranges on each chromosome.

PARAMETER	DESCRIPTION
ranges	Input genomic ranges. TYPE: Sequence[GenomicPosition]
genome	Optional chromosome lengths for boundary gaps. TYPE: dict[str, int] | None DEFAULT: None

RETURNS	DESCRIPTION
list[GenomicPosition]	Ranges representing uncovered regions.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def gaps(
    self, ranges: Sequence[GenomicPosition], *, genome: dict[str, int] | None = None
) -> list[GenomicPosition]:
    """Return gaps between ranges on each chromosome.

    Args:
        ranges: Input genomic ranges.
        genome: Optional chromosome lengths for boundary gaps.

    Returns:
        Ranges representing uncovered regions.
    """
    gr = self._gr(ranges)
    kwargs = {"seqlengths": genome} if genome else {}
    return self._pos(gr.gaps(**kwargs))

disjoin

disjoin(ranges)

Split ranges into non-overlapping disjoint intervals.

PARAMETER	DESCRIPTION
ranges	Input genomic ranges, possibly overlapping. TYPE: Sequence[GenomicPosition]

RETURNS	DESCRIPTION
list[GenomicPosition]	Maximal set of non-overlapping ranges.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def disjoin(self, ranges: Sequence[GenomicPosition]) -> list[GenomicPosition]:
    """Split ranges into non-overlapping disjoint intervals.

    Args:
        ranges: Input genomic ranges, possibly overlapping.

    Returns:
        Maximal set of non-overlapping ranges.
    """
    return self._pos(self._gr(ranges).disjoin())

nearest

nearest(query, subject, *, ignore_strand=True)

Find the nearest subject range for each query range.

PARAMETER	DESCRIPTION
query	Ranges to search for neighbors. TYPE: Sequence[GenomicPosition]
subject	Ranges to search within. TYPE: Sequence[GenomicPosition]
ignore_strand	If False, only consider same-strand ranges. TYPE: bool DEFAULT: True

RETURNS	DESCRIPTION
list[int | None]	List of indices into subject. None if no neighbor found.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def nearest(
    self,
    query: Sequence[GenomicPosition],
    subject: Sequence[GenomicPosition],
    *,
    ignore_strand: bool = True,
) -> list[int | None]:
    """Find the nearest subject range for each query range.

    Args:
        query: Ranges to search for neighbors.
        subject: Ranges to search within.
        ignore_strand: If False, only consider same-strand ranges.

    Returns:
        List of indices into subject. None if no neighbor found.
    """
    result = self._gr(query).nearest(self._gr(subject), ignore_strand=ignore_strand)
    # BiocPy returns -1 or None for no match; normalize to None
    return [None if (x is None or x < 0) else int(x) for x in result]

precede

precede(query, subject, *, ignore_strand=True)

Find the subject range immediately upstream of each query range.

PARAMETER	DESCRIPTION
query	Ranges to search for upstream neighbors. TYPE: Sequence[GenomicPosition]
subject	Ranges to search within. TYPE: Sequence[GenomicPosition]
ignore_strand	If False, only consider same-strand ranges. TYPE: bool DEFAULT: True

RETURNS	DESCRIPTION
list[int | None]	List of indices into subject. None if no upstream range found.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def precede(
    self,
    query: Sequence[GenomicPosition],
    subject: Sequence[GenomicPosition],
    *,
    ignore_strand: bool = True,
) -> list[int | None]:
    """Find the subject range immediately upstream of each query range.

    Args:
        query: Ranges to search for upstream neighbors.
        subject: Ranges to search within.
        ignore_strand: If False, only consider same-strand ranges.

    Returns:
        List of indices into subject. None if no upstream range found.
    """
    result = self._gr(query).precede(self._gr(subject), ignore_strand=ignore_strand)
    return [None if (x is None or x < 0) else int(x) for x in result]

follow

follow(query, subject, *, ignore_strand=True)

Find the subject range immediately downstream of each query range.

PARAMETER	DESCRIPTION
query	Ranges to search for downstream neighbors. TYPE: Sequence[GenomicPosition]
subject	Ranges to search within. TYPE: Sequence[GenomicPosition]
ignore_strand	If False, only consider same-strand ranges. TYPE: bool DEFAULT: True

RETURNS	DESCRIPTION
list[int | None]	List of indices into subject. None if no downstream range found.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def follow(
    self,
    query: Sequence[GenomicPosition],
    subject: Sequence[GenomicPosition],
    *,
    ignore_strand: bool = True,
) -> list[int | None]:
    """Find the subject range immediately downstream of each query range.

    Args:
        query: Ranges to search for downstream neighbors.
        subject: Ranges to search within.
        ignore_strand: If False, only consider same-strand ranges.

    Returns:
        List of indices into subject. None if no downstream range found.
    """
    result = self._gr(query).follow(self._gr(subject), ignore_strand=ignore_strand)
    return [None if (x is None or x < 0) else int(x) for x in result]

flank

flank(ranges, width, *, start=True, both=False)

Return flanking regions for each input range.

PARAMETER	DESCRIPTION
ranges	Input genomic ranges. TYPE: Sequence[GenomicPosition]
width	Width of the flanking region in bases. TYPE: int
start	If True, flank the start; if False, flank the end. TYPE: bool DEFAULT: True
both	If True, return flanks on both sides. TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
list[GenomicPosition]	Flanking ranges.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def flank(
    self,
    ranges: Sequence[GenomicPosition],
    width: int,
    *,
    start: bool = True,
    both: bool = False,
) -> list[GenomicPosition]:
    """Return flanking regions for each input range.

    Args:
        ranges: Input genomic ranges.
        width: Width of the flanking region in bases.
        start: If True, flank the start; if False, flank the end.
        both: If True, return flanks on both sides.

    Returns:
        Flanking ranges.
    """
    return self._pos(self._gr(ranges).flank(width=width, start=start, both=both))

resize

resize(ranges, width, *, fix='start')

Resize ranges to a fixed width anchored at start, end, or center.

PARAMETER	DESCRIPTION
ranges	Input genomic ranges. TYPE: Sequence[GenomicPosition]
width	Target width in bases. TYPE: int
fix	Anchor point — one of 'start', 'end', or 'center'. TYPE: str DEFAULT: 'start'

RETURNS	DESCRIPTION
list[GenomicPosition]	Resized ranges.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def resize(
    self, ranges: Sequence[GenomicPosition], width: int, *, fix: str = "start"
) -> list[GenomicPosition]:
    """Resize ranges to a fixed width anchored at start, end, or center.

    Args:
        ranges: Input genomic ranges.
        width: Target width in bases.
        fix: Anchor point — one of 'start', 'end', or 'center'.

    Returns:
        Resized ranges.
    """
    return self._pos(self._gr(ranges).resize(width=width, fix=fix))

shift

shift(ranges, shift)

Shift all ranges by a fixed number of bases.

PARAMETER	DESCRIPTION
ranges	Input genomic ranges. TYPE: Sequence[GenomicPosition]
shift	Bases to shift. Positive = downstream, negative = upstream. TYPE: int

RETURNS	DESCRIPTION
list[GenomicPosition]	Shifted ranges.

Source code in src/iris/adapters/driven/genomic_analysis/biocpy.py

def shift(self, ranges: Sequence[GenomicPosition], shift: int) -> list[GenomicPosition]:
    """Shift all ranges by a fixed number of bases.

    Args:
        ranges: Input genomic ranges.
        shift: Bases to shift. Positive = downstream, negative = upstream.

    Returns:
        Shifted ranges.
    """
    return self._pos(self._gr(ranges).shift(shift=shift))

pyhpo operations

iris.adapters.driven.genomic_analysis.pyhpo

PyHPO adapter for the HpoOperations port.

Bridges between iris Phenotype domain objects and the pyhpo library. Term resolution follows a two-step strategy: prefer hpo_id when available, fall back to name otherwise.

pyhpo initialises a module-level singleton on first Ontology() call. Subsequent PyHpoOperations instances reuse the same ontology in memory.

PyHpoOperations

PyHpoOperations(*, data_directory=None)

Bases: HpoOperations

HpoOperations adapter backed by the pyhpo library.

On construction the HPO ontology is loaded into memory (once per process — pyhpo is a singleton). Subsequent instantiations reuse the already-loaded data without re-reading from disk.

PARAMETER	DESCRIPTION
data_directory	Path to a directory containing HPO data files. When None (default), pyhpo uses its own bundled data. TYPE: str | None DEFAULT: None

Example

from iris.adapters.driven.genomic_analysis.pyhpo import PyHpoOperations
from iris.domain.genomics.services import HpoService

service = HpoService(ops=PyHpoOperations())
hits = service.search("Seizure", limit=5)
score = service.similarity(hits[0], hits[1])

Source code in src/iris/adapters/driven/genomic_analysis/pyhpo.py

def __init__(self, *, data_directory: str | None = None) -> None:
    onto, bundled = _require_pyhpo()
    onto(data_directory or bundled)
    self._onto = onto

PyHpoEnrichmentModel

PyHpoEnrichmentModel(*, data_directory=None)

Bases: PhenotypeEnrichmentModel

PhenotypeEnrichmentModel adapter backed by pyhpo's EnrichmentModel.

On construction the HPO ontology is loaded into memory (once per process — pyhpo is a singleton). The 'gene' and 'omim' EnrichmentModel instances are stateful and built lazily on first use, then cached for the lifetime of the adapter.

PARAMETER	DESCRIPTION
data_directory	Path to a directory containing HPO data files. When None (default), pyhpo uses its own bundled data. TYPE: str | None DEFAULT: None

Example

from iris.adapters.driven.genomic_analysis.pyhpo import PyHpoEnrichmentModel

model = PyHpoEnrichmentModel()
hits = model.gene_enrichment(patient_phenotypes, limit=20)
top_symbols = [h.gene.symbol for h in hits if h.is_significant()]

Source code in src/iris/adapters/driven/genomic_analysis/pyhpo.py

def __init__(self, *, data_directory: str | None = None) -> None:
    onto, bundled = _require_pyhpo()
    onto(data_directory or bundled)
    self._onto = onto
    self._gene_model: Any = None
    self._disease_model: Any = None

Genotype processors

iris.adapters.driven.genomic_analysis.processors

Concrete GenotypeProcessor implementations for cohort-level analyses.

AncestryStratifiedAFProcessor

AncestryStratifiedAFProcessor(
    sample_amr,
    *,
    valid_mask=None,
    n_deciles=10,
    p_threshold=0.05,
    min_ac_per_decile=5,
    min_valid_deciles=3,
)

Compute per-decile allele frequencies and an ancestry linear regression.

For each variant computes:

• af_decile: list of length n_deciles with AF per AMR quantile bin.
• amr_beta: regression slope AF ~ AMR decile midpoint.
• amr_se: standard error of the slope.
• amr_pvalue: two-sided p-value of the slope.
• amr_r2: coefficient of determination R².

PARAMETER	DESCRIPTION
sample_amr	float array (n_samples,) with AMR ancestry proportion per sample, aligned to the Zarr sample_id dimension. TYPE: Any
n_deciles	number of quantile bins. Default 10. TYPE: int DEFAULT: 10
p_threshold	significance threshold for downstream filtering. Default 0.05. TYPE: float DEFAULT: 0.05
min_ac_per_decile	minimum AC in a decile to include it in the regression. Default 5. TYPE: int DEFAULT: 5
min_valid_deciles	minimum non-NaN deciles required to run the regression. Default 3. TYPE: int DEFAULT: 3

Example

import numpy as np
import pandas as pd

covars = pd.read_csv("ancestry.csv", index_col="sample_id")
amr = covars.loc[sample_order, "AMR"].values

processor = AncestryStratifiedAFProcessor(sample_amr=amr)

with ZarrVariantStore(path=ZARR_PATH, processor=processor) as store:
    variants = store.fetch(region)

for va in variants:
    print(va.metadata["af_decile"])  # list[float | None]
    print(va.metadata["amr_pvalue"])  # float

Initialize the processor.

PARAMETER	DESCRIPTION
sample_amr	float array (n_valid_samples,) with AMR ancestry proportion, ALREADY FILTERED to only the samples that have a known ancestry value (same length and order as the True entries in valid_mask, if provided). TYPE: Any
valid_mask	optional bool array (n_zarr_samples,) — True for samples present in the Zarr store that also have a known ancestry value. Samples with False are excluded from every computation (AF, deciles, regression). When None, all Zarr samples are assumed valid and sample_amr must cover all of them. TYPE: Any | None DEFAULT: None
n_deciles	number of quantile bins. Default 10. TYPE: int DEFAULT: 10
p_threshold	significance threshold for downstream filtering. TYPE: float DEFAULT: 0.05
min_ac_per_decile	minimum AC in a decile to include it. TYPE: int DEFAULT: 5
min_valid_deciles	minimum non-NaN deciles to run regression. TYPE: int DEFAULT: 3

Source code in src/iris/adapters/driven/genomic_analysis/processors.py

def __init__(
    self,
    sample_amr: Any,
    *,
    valid_mask: Any | None = None,
    n_deciles: int = 10,
    p_threshold: float = 0.05,
    min_ac_per_decile: int = 5,
    min_valid_deciles: int = 3,
) -> None:
    """Initialize the processor.

    Args:
        sample_amr: float array (n_valid_samples,) with AMR ancestry
            proportion, ALREADY FILTERED to only the samples that have
            a known ancestry value (same length and order as the True
            entries in valid_mask, if provided).
        valid_mask: optional bool array (n_zarr_samples,) — True for
            samples present in the Zarr store that also have a known
            ancestry value. Samples with False are excluded from every
            computation (AF, deciles, regression). When None, all Zarr
            samples are assumed valid and sample_amr must cover all of
            them.
        n_deciles: number of quantile bins. Default 10.
        p_threshold: significance threshold for downstream filtering.
        min_ac_per_decile: minimum AC in a decile to include it.
        min_valid_deciles: minimum non-NaN deciles to run regression.
    """
    self._amr = np.asarray(sample_amr, dtype=np.float64)
    self._valid_mask = np.asarray(valid_mask, dtype=bool) if valid_mask is not None else None
    self._n_deciles = n_deciles
    self._p_threshold = p_threshold
    self._min_ac = min_ac_per_decile
    self._min_valid = min_valid_deciles

    quantiles = np.linspace(0, 100, n_deciles + 1)
    breakpoints = np.percentile(self._amr, quantiles)
    self._decile_idx: Any = np.clip(
        np.digitize(self._amr, bins=breakpoints[1:-1]), 0, n_deciles - 1
    )
    self._decile_midpoints: Any = np.array(
        [
            (
                self._amr[self._decile_idx == d].mean()
                if (self._decile_idx == d).any()
                else np.nan
            )
            for d in range(n_deciles)
        ],
        dtype=np.float64,
    )

HomozygousByQuintileProcessor

HomozygousByQuintileProcessor(
    sample_quintile, *, valid_mask=None, n_quintiles=5
)

Bases: GenotypeProcessor

Cuenta homocigotos alt (2 alelos alternativos en la misma muestra) por quintil.

Deprecated: usar AltCountProcessor(sample_strata=[SampleStratum(f"q{i}", mask), ...]) — mismo resultado (hom_qX/called_qX equivalen a all__qX__n_hom/ all__qX__n_called) más estratificación por variante y acumulación por muestra (summary_stats()).

Para cada variante, agrupa las muestras por quintil (1..n_quintiles — p.ej. quintil de ancestría AMR) y cuenta cuántas son homocigotas para el alelo alternativo (ambos alelos == alt_idx) frente a cuántas tienen genotipo llamado en ese quintil.

Agrega al metadata de cada variante

hom_q{i} — homocigotos alt en el quintil i called_q{i} — muestras con genotipo llamado (no missing) en el quintil i hom_rate_q{i} — hom_q{i} / called_q{i} (NaN si called_q{i} == 0) hom_total — homocigotos alt sumados sobre todos los quintiles válidos called_total — llamados sumados sobre todos los quintiles válidos hom_rate_total — hom_total / called_total (NaN si called_total == 0)

PARAMETER	DESCRIPTION
sample_quintile	int array (n_valid_samples,) con el quintil (1.. n_quintiles) de cada muestra VÁLIDA, ya alineado y filtrado al mismo orden que las entradas True de valid_mask — mismo contrato que AncestryStratifiedAFProcessor.sample_amr. TYPE: Any
valid_mask	bool array opcional (n_zarr_samples,) — True para las muestras del store que tienen un quintil conocido. Las muestras con False quedan fuera de todo el cómputo (nunca cuentan como homocigotas ni como llamadas). Si es None, se asume que todas las muestras del store son válidas y sample_quintile debe cubrirlas todas, en el mismo orden. TYPE: Any | None DEFAULT: None
n_quintiles	número de quintiles. Default 5. TYPE: int DEFAULT: 5

Example

import csv
import numpy as np

quintile_by_sample: dict[str, int] = {}
with open("mcps_ancestry_quintiles.csv") as fh:
    for row in csv.DictReader(fh):
        q = (row["amr_quintil"] or "").strip()
        if q.isdigit():
            quintile_by_sample[row["sample_code"]] = int(q)

with ZarrVariantStore(path=ZARR_PATH) as store:
    valid, quintiles = [], []
    for sample in store.sample_ids:
        q = quintile_by_sample.get(sample)
        valid.append(q is not None)
        if q is not None:
            quintiles.append(q)

    sq, vm = np.array(quintiles), np.array(valid)
    store.processor = HomozygousByQuintileProcessor(sample_quintile=sq, valid_mask=vm)
    variants = store.fetch(region)

for va in variants:
    print(va.metadata["hom_q1"], va.metadata["hom_rate_q1"])

Initialize the processor.

PARAMETER	DESCRIPTION
sample_quintile	int array (n_valid_samples,), ver docstring de clase. TYPE: Any
valid_mask	bool array opcional (n_zarr_samples,), ver docstring de clase. TYPE: Any | None DEFAULT: None
n_quintiles	número de quintiles. Default 5. TYPE: int DEFAULT: 5

Source code in src/iris/adapters/driven/genomic_analysis/processors.py

def __init__(
    self, sample_quintile: Any, *, valid_mask: Any | None = None, n_quintiles: int = 5
) -> None:
    """Initialize the processor.

    Args:
        sample_quintile: int array (n_valid_samples,), ver docstring de clase.
        valid_mask: bool array opcional (n_zarr_samples,), ver docstring de clase.
        n_quintiles: número de quintiles. Default 5.
    """
    warnings.warn(
        "HomozygousByQuintileProcessor está deprecado — usar AltCountProcessor "
        "con sample_strata (ver docstring de clase).",
        DeprecationWarning,
        stacklevel=2,
    )
    self._quintile = np.asarray(sample_quintile, dtype=np.int64)
    self._valid_mask = np.asarray(valid_mask, dtype=bool) if valid_mask is not None else None
    self._n_quintiles = n_quintiles

PerIndividualAltCountProcessor

PerIndividualAltCountProcessor()

Bases: GenotypeProcessor

Calcula la distribución de alelos alt por individuo sobre un bloque.

Deprecated: usar AltCountProcessor() (sin strata) — mismo resultado (all__all__alt_counts_median/all__all__vars_median, etc.) más estratificación opcional por muestra/variante y acumulación por muestra (summary_stats()).

Agrega al metadata de cada variante

alt_counts_median — mediana de (alelos alt por individuo) en la cohorte alt_counts_iqr — IQR (p75 - p25) vars_median — mediana de variantes portadas por individuo vars_iqr — IQR

Nota: 'alt_counts' = 0/½ por individuo; 'vars' = 1 si portador (het o hom), 0 si no.

Initialize the processor.

Source code in src/iris/adapters/driven/genomic_analysis/processors.py

def __init__(self) -> None:
    """Initialize the processor."""
    warnings.warn(
        "PerIndividualAltCountProcessor está deprecado — usar AltCountProcessor() "
        "sin strata (ver docstring de clase).",
        DeprecationWarning,
        stacklevel=2,
    )

SampleStratum dataclass

SampleStratum(name, mask)

Un grupo de muestras sobre el que estratificar (p.ej. un quintil de ancestría).

ATTRIBUTE	DESCRIPTION
name	nombre del stratum — usado en el prefijo de las keys de salida. TYPE: str
mask	bool array (n_samples_zarr,) — True para las muestras del store que pertenecen a este stratum, en el mismo orden que la dimensión de muestras del Zarr (ver ZarrVariantStore.sample_ids). TYPE: ndarray

VariantStratum dataclass

VariantStratum(name, cpra_set)

Un grupo de variantes sobre el que estratificar (p.ej. un panel de genes).

ATTRIBUTE	DESCRIPTION
name	nombre del stratum — usado en el prefijo de las keys de salida. TYPE: str
cpra_set	set de identificadores CHR:POS:REF:ALT que pertenecen a este stratum, o None para incluir todas las variantes del bloque sin filtrar. TYPE: set[str] | None

AltCountProcessor

AltCountProcessor(
    sample_strata=None,
    variant_strata=None,
    *,
    compute_percentiles=True,
)

Bases: GenotypeProcessor

Conteo de alelos alt estratificado por muestra y por variante, ortogonal.

Generaliza a PerIndividualAltCountProcessor (sin estratificar) y a HomozygousByQuintileProcessor (estratificado por muestra) bajo una sola implementación con dos ejes de estratificación independientes:

• sample_strata: agrupa las muestras (p.ej. quintiles de ancestría). Por defecto (None) usa un único stratum "all" con todas las muestras del bloque.
• variant_strata: agrupa las variantes por identidad CPRA vía context["cpra"] (p.ej. un panel de genes). Por defecto (None) usa un único stratum "all" sin filtrar.

Para cada combinación (variant_stratum, sample_stratum) agrega al metadata de cada variante, bajo el prefijo {vs_name}__{ss_name}__: n_het, n_hom, n_called, af, alt_counts_median, alt_counts_q25, alt_counts_q75, vars_median, vars_q25, vars_q75 — todas NaN en las variantes que no pertenecen al variant_stratum. Las últimas 6 (percentiles/medianas) se pueden omitir con compute_percentiles=False — son las más caras de calcular y, con muchas combinaciones de strata, las que más pesan en el metadata por variante cuando el caller no las necesita (ver __init__).

Además acumula, por bloque, métricas por muestra recuperables al final vía :meth:summary_stats:

• n_variants: variantes portadas (het u hom) por la muestra, sumadas sobre todos los bloques — no depende de la agrupación por gen.
• n_genes_affected: 1 por gen si la muestra porta >= 1 variante del variant_stratum en ese gen (sumado sobre todos los genes vistos).
• n_genes_biallelic: 1 por gen si la muestra es hom en alguna variante del gen, o portadora de >= 2 variantes het distintas del gen (en cualquier combinación de bloques de ese gen).

Identidad de "gen": context["gene"] — que ZarrVariantStore completa automáticamente con region.name (la columna 4 del BED que generó esa región) si el caller no lo puso ya. Cuando el BED tiene una fila por gen (*.genes.bed), cada bloque ES un gen y el resultado es idéntico al de antes. Cuando el BED tiene varias filas por gen (*.exons.bed — un exón por fila), los bloques con el mismo context["gene"] se combinan: portador en dos exones distintos del mismo gen sigue contando n_genes_affected una sola vez, y compound-het entre exones (variante het en el exón 2, otra variante het distinta en el exón 7) sí se detecta como bialélico, sin volver a incluir los intrones entre exones en el cómputo (cada bloque sigue siendo solo la región exacta que le pasó el caller).

Si context["gene"] es None (BED sin columna de nombre, o el caller no la puso), cada bloque vuelve a tratarse como su propio "gen" — comportamiento histórico, sin agrupación.

n_genes_affected/n_genes_biallelic en _acc/summary_stats() se derivan del estado por gen recién al llamar :meth:finalize (que :meth:summary_stats invoca automáticamente) — no durante el escaneo, porque los bloques de un mismo gen pueden llegar en cualquier orden o en hilos distintos. Si leés processor._acc directamente en vez de summary_stats() (como hace 04_gene_stats.py), llamá processor.finalize() primero.

PARAMETER	DESCRIPTION
sample_strata	lista de SampleStratum, o None para un único stratum "all" con todas las muestras del bloque. TYPE: list[SampleStratum] | None DEFAULT: None
variant_strata	lista de VariantStratum, o None para un único stratum "all" sin filtrar. TYPE: list[VariantStratum] | None DEFAULT: None

Example

strata = []
for q in range(1, 6):
    strata.append(SampleStratum(f"q{q}", sample_quintile == q))

processor = AltCountProcessor(sample_strata=strata)
with ZarrVariantStore(path=ZARR_PATH, processor=processor) as store:
    variants = store.fetch(region)

for va in variants:
    print(va.metadata["all__q1__n_hom"])

stats = processor.summary_stats()
print(stats["all"]["q1"]["n_genes_affected"]["median"])

Initialize the processor.

PARAMETER	DESCRIPTION
sample_strata	ver docstring de clase. TYPE: list[SampleStratum] | None DEFAULT: None
variant_strata	ver docstring de clase. TYPE: list[VariantStratum] | None DEFAULT: None
compute_percentiles	si es False, omite alt_counts_median/q25/q75 y vars_median/q25/q75 del dict que retorna __call__ (queda solo n_het, n_hom, n_called, af). No afecta :meth:summary_stats — la acumulación por muestra (_acc) es independiente de estas métricas por variante. Poné esto en False cuando el caller no lee esas 6 columnas (p.ej. solo le interesa el acumulador, o solo un subconjunto de las métricas) — evita 3 llamadas a np.nanpercentile por combinación (variant_stratum, sample_stratum) y bloque, que son las más caras de las 10 métricas. TYPE: bool DEFAULT: True

Source code in src/iris/adapters/driven/genomic_analysis/processors.py

def __init__(
    self,
    sample_strata: list[SampleStratum] | None = None,
    variant_strata: list[VariantStratum] | None = None,
    *,
    compute_percentiles: bool = True,
) -> None:
    """Initialize the processor.

    Args:
        sample_strata: ver docstring de clase.
        variant_strata: ver docstring de clase.
        compute_percentiles: si es ``False``, omite ``alt_counts_median/q25/q75``
            y ``vars_median/q25/q75`` del dict que retorna ``__call__`` (queda
            solo ``n_het``, ``n_hom``, ``n_called``, ``af``). No afecta
            :meth:`summary_stats` — la acumulación por muestra (``_acc``) es
            independiente de estas métricas por variante. Poné esto en ``False``
            cuando el caller no lee esas 6 columnas (p.ej. solo le interesa el
            acumulador, o solo un subconjunto de las métricas) — evita 3
            llamadas a ``np.nanpercentile`` por combinación
            ``(variant_stratum, sample_stratum)`` y bloque, que son las más
            caras de las 10 métricas.
    """
    self._sample_strata = sample_strata
    self._variant_strata: list[VariantStratum] = (
        variant_strata if variant_strata is not None else [VariantStratum("all", None)]
    )
    self._sample_mask_by_name: dict[str, np.ndarray] = {}
    self._acc: dict[str, dict[str, dict[str, np.ndarray]]] = {}
    self._compute_percentiles = compute_percentiles
    # Estado crudo por gen — [vs_name][ss_name][gene] -> {"carrier",
    # "hom", "het_count"} — combinado entre todos los bloques (exones)
    # que compartan el mismo context["gene"]. Ver finalize().
    self._gene_state: dict[str, dict[str, dict[str, dict[str, np.ndarray]]]] = {}
    self._gene_lock = threading.Lock()

finalize

finalize()

Deriva n_genes_affected/n_genes_biallelic finales desde el estado por gen.

Solo hace falta llamarlo explícitamente si leés processor._acc directamente (p.ej. 04_gene_stats.py) — summary_stats() ya lo invoca. Idempotente: siempre recalcula desde self._gene_state en vez de incrementar, así que llamarlo varias veces no duplica conteos. No toca las combinaciones (vs, ss) que nunca recibieron un context["gene"] — esas ya quedaron acumuladas de inmediato en _accumulate().

Source code in src/iris/adapters/driven/genomic_analysis/processors.py

def finalize(self) -> None:
    """Deriva ``n_genes_affected``/``n_genes_biallelic`` finales desde el estado por gen.

    Solo hace falta llamarlo explícitamente si leés ``processor._acc``
    directamente (p.ej. ``04_gene_stats.py``) — ``summary_stats()`` ya lo
    invoca. Idempotente: siempre recalcula desde ``self._gene_state`` en
    vez de incrementar, así que llamarlo varias veces no duplica conteos.
    No toca las combinaciones (vs, ss) que nunca recibieron un
    ``context["gene"]`` — esas ya quedaron acumuladas de inmediato en
    ``_accumulate()``.
    """
    with self._gene_lock:
        for vs_name, by_ss in self._gene_state.items():
            for ss_name, genes in by_ss.items():
                acc = self._acc.get(vs_name, {}).get(ss_name)
                if acc is None or not genes:
                    continue
                n_samples = len(acc["n_genes_affected"])
                affected_total = np.zeros(n_samples, dtype=np.int64)
                biallelic_total = np.zeros(n_samples, dtype=np.int64)
                for state in genes.values():
                    affected_total += state["carrier"].astype(np.int64)
                    biallelic_total += (state["hom"] | (state["het_count"] >= 2)).astype(
                        np.int64
                    )
                acc["n_genes_affected"] = affected_total
                acc["n_genes_biallelic"] = biallelic_total

summary_stats

summary_stats()

Distribución (median/q25/q75/mean/max) de las métricas por muestra.

Llama :meth:finalize primero — ver su docstring y el de la clase sobre agrupación por gen.

Solo considera, para cada SampleStratum, las muestras que pertenecen a ese stratum (ss.mask) — nunca las que quedaron fuera de todo stratum de muestra.

RETURNS	DESCRIPTION
dict[str, dict[str, dict[str, dict[str, float]]]]	{vs_name: {ss_name: {metric: {stat: value}}}}.

Source code in src/iris/adapters/driven/genomic_analysis/processors.py

def summary_stats(self) -> dict[str, dict[str, dict[str, dict[str, float]]]]:
    """Distribución (median/q25/q75/mean/max) de las métricas por muestra.

    Llama :meth:`finalize` primero — ver su docstring y el de la clase
    sobre agrupación por gen.

    Solo considera, para cada ``SampleStratum``, las muestras que
    pertenecen a ese stratum (``ss.mask``) — nunca las que quedaron
    fuera de todo stratum de muestra.

    Returns:
        ``{vs_name: {ss_name: {metric: {stat: value}}}}``.
    """
    self.finalize()
    out: dict[str, dict[str, dict[str, dict[str, float]]]] = {}
    for vs_name, by_ss in self._acc.items():
        out[vs_name] = {}
        for ss_name, metrics in by_ss.items():
            mask = self._sample_mask_by_name.get(ss_name)
            valid_metrics: dict[str, dict[str, float]] = {}
            for metric_name, arr in metrics.items():
                valid = arr[mask] if mask is not None else arr
                valid_metrics[metric_name] = _distribution(valid)
            out[vs_name][ss_name] = valid_metrics
    return out