Stores & Pipeline

Adapters for cohort-scale variant extraction and parallel annotation pipelines.

ZarrVariantStore

Cohort store backed by a vcf2zarr Zarr v3 archive. Implements the VariantAnnotationStore port.

Key design points

• variant_position values are 1-based (VCF native, as written by vcf2zarr without --zero-based).
• open() loads the full variant_contig array once and builds per-contig start/end offsets via numpy.searchsorted, then discards the array.
• Per-chromosome position vectors are pre-loaded into _chrom_pos_cache for O(1) mask computation inside fetch().
• fetch() is thread-safe: it performs only read operations on pre-loaded numpy arrays and issues independent fancy-index reads against the Zarr store (Zarr's LocalStore supports concurrent reads).

Usage

from iris.adapters.driven.stores import ZarrVariantStore

with ZarrVariantStore(path="/data/wes_cohort.zarr") as store:
    variants = store.fetch(region)

Genotype processors

ZarrVariantStore accepts an optional processor (a GenotypeProcessor implementation) and processor_context dict. When set, the processor is called once per fetch block per ALT allele with the raw genotype arrays (n_variants, n_samples, ploidy), and its returned dict is merged into each VariantAnnotation.metadata.data. Built-in genotype counts (n_hom, n_het, …) are always computed regardless of whether a processor is supplied. See Genomic analysis for concrete implementations (AncestryStratifiedAFProcessor, PerIndividualAltCountProcessor).

iris.adapters.driven.stores.zarr

ZarrVariantStore — vcf2zarr cohort store adapter.

Supports user-defined GenotypeProcessor callables injected at construction time, enabling cohort-level statistics, stratified allele frequencies, and regression models without subclassing.

ZarrVariantStore

Bases: VariantAnnotationStore

Cohort-scale variant store backed by a vcf2zarr Zarr v3 archive.

Designed for variant_position values that are 1-based (VCF native, as written by vcf2zarr without the --zero-based flag).

The store must be opened before use — either explicitly via :meth:open or through the context-manager protocol::

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

After open():

• Contig offsets are built once from the full variant_contig array using numpy.searchsorted.
• Per-chromosome position arrays are pre-loaded and cached for O(1) mask computation in :meth:fetch.

:meth:fetch is thread-safe: it performs only read operations on the pre-loaded numpy arrays and issues independent fancy-index reads against the underlying Zarr store (Zarr's LocalStore supports concurrent reads).

ATTRIBUTE	DESCRIPTION
path	Path to the Zarr store directory or .zarr file. TYPE: Path
assembly	Genome assembly label stored in created :class:~iris.domain.genomics.objects.GenomicPosition objects. TYPE: str
source	Population / source label used in the generated :class:~iris.domain.genomics.objects.AlleleFrequency objects. TYPE: str

processor class-attribute instance-attribute

processor = None

Optional user-defined processor called for each variant block.

When set, the processor receives the raw genotype arrays (n_variants, n_samples, ploidy) and returns a dict of computed values that are merged into VariantAnnotation.metadata.data. The built-in genotype counts (n_hom, n_het, etc.) are always computed regardless of whether a processor is supplied.

See :class:~iris.domain.genomics.ports.driven.genotype_processor.GenotypeProcessor.

processor_context class-attribute instance-attribute

processor_context = field(factory=dict)

Arbitrary context dict forwarded to the processor on every call.

Use this to pass sample covariates, phenotypes, or configuration that the processor needs but that are not part of the genotype arrays.

open

open()

Open the Zarr store and build in-memory indices.

Loads the full variant_contig array (~35 MB for a 79k cohort) once to build per-contig start / end offsets, then pre-loads each per-chromosome position array into _chrom_pos_cache. The large variant_contig array is discarded after the index is built.

RETURNS	DESCRIPTION
ZarrVariantStore	Self, so this can be used fluently: store.open().

RAISES	DESCRIPTION
ImportError	If zarr or numpy are not installed.

Source code in src/iris/adapters/driven/stores/zarr.py

def open(self) -> ZarrVariantStore:
    """Open the Zarr store and build in-memory indices.

    Loads the full ``variant_contig`` array (~35 MB for a 79k cohort)
    once to build per-contig start / end offsets, then pre-loads each
    per-chromosome position array into ``_chrom_pos_cache``.  The large
    ``variant_contig`` array is discarded after the index is built.

    Returns:
        Self, so this can be used fluently: ``store.open()``.

    Raises:
        ImportError: If ``zarr`` or ``numpy`` are not installed.
    """
    try:
        import numpy as np  # pyright: ignore[reportMissingImports]
        import zarr  # pyright: ignore[reportMissingImports]
    except ImportError as exc:
        raise ImportError(
            "zarr and numpy are required: uv pip install 'iris[genomics]'"
        ) from exc

    self._store = zarr.open(str(self.path), mode="r")
    contig_id = self._store["contig_id"][:]
    self._contig_names = [
        HUMAN_CHROMOSOMES.resolve(str(c)) for c in contig_id if str(c) in HUMAN_CHROMOSOMES
    ]
    self._contig_name_to_idx = {c: i for i, c in enumerate(self._contig_names)}

    # Build start/end offsets for each contig index using searchsorted.
    v_contig_all = self._store["variant_contig"][:]
    n_contigs = len(self._contig_names)
    idx = np.arange(n_contigs)
    self._contig_start = np.searchsorted(v_contig_all, idx, side="left").astype(int)
    self._contig_end = np.searchsorted(v_contig_all, idx, side="right").astype(int)
    del v_contig_all

    # Pre-load per-chromosome position vectors.
    for chrom, ci in self._contig_name_to_idx.items():
        s, e = int(self._contig_start[ci]), int(self._contig_end[ci])
        if s >= e:
            continue
        pos = self._store["variant_position"][s:e]
        self._chrom_pos_cache[chrom] = (pos, s)

    self._n_samples = len(self._store["sample_id"][:]) if "sample_id" in self._store else 0
    return self

close

close()

Release in-memory caches and close the Zarr store handle.

Source code in src/iris/adapters/driven/stores/zarr.py

def close(self) -> None:
    """Release in-memory caches and close the Zarr store handle."""
    self._store = None
    self._chrom_pos_cache = {}
    self._contig_start = None
    self._contig_end = None

fetch

fetch(region)

Return all variants whose positions fall within region.

region.start and region.end must both be non-None and are interpreted as 1-based inclusive coordinates matching the variant_position encoding in the store.

The method opens the store lazily if :meth:open has not been called. It is safe to call concurrently from multiple threads.

PARAMETER	DESCRIPTION
region	1-based inclusive genomic window. TYPE: GenomicPosition

RETURNS	DESCRIPTION
list[VariantAnnotation]	Possibly-empty list of
list[VariantAnnotation]	class:~iris.domain.genomics.entities.VariantAnnotation.

RAISES	DESCRIPTION
ValueError	If region.start or region.end is None.
ImportError	If zarr or numpy are not installed.

Source code in src/iris/adapters/driven/stores/zarr.py

def fetch(self, region: GenomicPosition) -> list[VariantAnnotation]:
    """Return all variants whose positions fall within *region*.

    ``region.start`` and ``region.end`` must both be non-``None`` and are
    interpreted as **1-based inclusive** coordinates matching the
    ``variant_position`` encoding in the store.

    The method opens the store lazily if :meth:`open` has not been called.
    It is safe to call concurrently from multiple threads.

    Args:
        region: 1-based inclusive genomic window.

    Returns:
        Possibly-empty list of
        :class:`~iris.domain.genomics.entities.VariantAnnotation`.

    Raises:
        ValueError: If ``region.start`` or ``region.end`` is ``None``.
        ImportError: If ``zarr`` or ``numpy`` are not installed.
    """
    if region.start is None or region.end is None:
        raise ValueError("region.start and region.end must be non-None for ZarrVariantStore")

    if self._store is None:
        self.open()

    try:
        import numpy as np  # pyright: ignore[reportMissingImports]
    except ImportError as exc:
        raise ImportError("numpy is required: uv pip install 'iris[genomics]'") from exc

    chrom = HUMAN_CHROMOSOMES.resolve(region.chromosome)
    if chrom not in self._chrom_pos_cache:
        return []

    chrom_pos, global_start = self._chrom_pos_cache[chrom]

    # 1-based inclusive position filter.
    mask = (chrom_pos >= region.start) & (chrom_pos <= region.end)
    local_indices: Any = np.where(mask)[0]
    if len(local_indices) == 0:
        return []

    global_indices = local_indices + global_start
    positions = chrom_pos[local_indices]

    # Load alleles and genotypes only for the filtered variants.
    v_allele = self._store["variant_allele"][global_indices]
    v_id = self._store["variant_id"][global_indices] if "variant_id" in self._store else None

    gt = self._store["call_genotype"][global_indices]  # (n_filt, n_samples, 2)
    gt_mask = (
        self._store["call_genotype_mask"][global_indices]
        if "call_genotype_mask" in self._store
        else None
    )

    # Run the user processor once per block (all variants, all samples)
    # for each ALT allele index present in this block.
    # Results are keyed by alt_idx so we can look them up per variant.
    # processor_block_results[alt_idx][field] -> array(n_filt,)
    processor_block_results: dict[int, dict[str, Any]] = {}
    if self.processor is not None:
        # Determine all alt indices present in this block
        max_alts = v_allele.shape[1] - 1  # allele columns minus ref
        for aidx in range(1, max_alts + 1):
            proc_out = self.processor(
                gt, gt_mask, alt_idx=aidx, context=self.processor_context
            )
            processor_block_results[aidx] = proc_out

    results: list[VariantAnnotation] = []

    for local_i in range(len(local_indices)):
        alleles_row = v_allele[local_i]
        ref = str(alleles_row[0])
        alts = [str(a) for a in alleles_row[1:] if str(a) not in ("", ".")]
        if not alts:
            continue

        pos = int(positions[local_i])  # already 1-based

        vid_raw = str(v_id[local_i]) if v_id is not None else None
        if vid_raw in ("", "."):
            vid_raw = None

        # Slice genotype arrays for this one variant.
        gt_v = gt[local_i : local_i + 1]  # shape (1, n_samples, 2)
        gt_mask_v = gt_mask[local_i : local_i + 1] if gt_mask is not None else None

        for alt_idx, alt in enumerate(alts, start=1):
            n_hom_ref_arr, n_het_arr, n_hom_alt_arr, n_missing_arr = (
                _genotype_counts_vectorized(gt_v, gt_mask_v, alt_idx)
            )
            gc = GenotypeCount(
                n_hom_ref=int(n_hom_ref_arr[0]),
                n_het=int(n_het_arr[0]),
                n_hom_alt=int(n_hom_alt_arr[0]),
                n_missing=int(n_missing_arr[0]),
                allele=alt,
            )
            af = gc.to_allele_frequency(self.source)

            vid = vid_raw or f"{chrom}-{pos}-{ref}-{alt}"
            dbsnp_id = vid_raw if vid_raw and vid_raw.startswith("rs") else None

            # Base metadata — always present
            meta: dict[str, Any] = {
                "n_hom": int(n_hom_alt_arr[0]),
                "n_het": int(n_het_arr[0]),
                "n_hom_ref": int(n_hom_ref_arr[0]),
                "n_missing": int(n_missing_arr[0]),
                "n_samples": self._n_samples,
                "ac": gc.ac,
                "an": gc.an,
                "af": (gc.ac / gc.an) if gc.an else None,
            }

            # Merge processor results for this variant and alt index
            if alt_idx in processor_block_results:
                for field, arr in processor_block_results[alt_idx].items():
                    if hasattr(arr, "__len__") and hasattr(arr, "__getitem__"):
                        # Array — extract the value for this variant
                        val = arr[local_i]
                        meta[field] = val.tolist() if hasattr(val, "tolist") else float(val)
                    else:
                        # Scalar — apply to all variants
                        meta[field] = arr

            results.append(
                VariantAnnotation(
                    variant=Variant(
                        variant_type=_infer_type(ref, alt),
                        nomenclature=VariantNomenclature.from_key(
                            (chrom, pos, ref, alt),
                            dbsnp_id=dbsnp_id,
                            extra_ids=[vid] if vid else None,
                        ),
                        alleles=Allele(reference=ref, alternate=(alt,)),
                        position=GenomicPosition(
                            chromosome=chrom,
                            assembly=self.assembly,
                            start=pos,
                            end=pos + len(ref) - 1,
                        ),
                    ),
                    allele_frequencies=(af,),
                    metadata=Metadata(source=self.source, data=meta),
                )
            )

    return results

VariantExtractionPipeline

Orchestrates parallel variant extraction followed by sequential annotation.

Two-phase execution

1. Extraction — the store's fetch() is called once per region using a ThreadPoolExecutor. Because Zarr stores release the GIL for array reads, multi-threading yields near-linear speedups on I/O-bound workloads.
2. Annotation — collected VariantAnnotation objects are passed sequentially through each annotator. Annotators run in serial to avoid saturating tabix file handles or in-memory lookup tables.

Usage

from iris.adapters.driven.pipeline import VariantExtractionPipeline
from iris.adapters.driven.stores import ZarrVariantStore
from iris.adapters.driven.annotators.clinvar import ClinvarAnnotator

annotators = (ClinvarAnnotator(path=CLINVAR_TSV),)
with ZarrVariantStore(path=ZARR_PATH) as store:
    pipeline = VariantExtractionPipeline(
        store, annotators, workers=28, show_progress=True
    )
    variants = pipeline.run(regions=regions_from_bed)

iris.adapters.driven.pipeline

VariantExtractionPipeline — parallel extraction and sequential annotation.

VariantExtractionPipeline

Orchestrates parallel variant extraction followed by sequential annotation.

The pipeline executes in two phases:

1. Extraction — :attr:store is queried once per region using a :class:~concurrent.futures.ThreadPoolExecutor with up to :attr:workers threads. Because Zarr stores release the GIL for array reads, multi-threading yields near-linear speedups.

2. Annotation — the collected :class:~iris.domain.genomics.entities.VariantAnnotation objects are passed through each annotator in :attr:annotators sequentially. Annotators run in serial to avoid saturating the I/O devices they read from (tabix files, in-memory lookup tables, etc.).

Example::

annotators = (ClinvarAnnotator(path=CLINVAR_TSV),)
with ZarrVariantStore(path=ZARR_PATH) as store:
    pipeline = VariantExtractionPipeline(store, annotators, workers=28)
    variants = pipeline.run(regions=regions_from_bed)

ATTRIBUTE	DESCRIPTION
store	Opened :class:~iris.domain.genomics.ports.driven.variant_store.VariantAnnotationStore to fetch variants from. TYPE: VariantAnnotationStore
annotators	Zero or more :class:~iris.domain.genomics.ports.driven.annotators.VariantAnnotator instances applied in order after extraction. TYPE: tuple[VariantAnnotator, ...]
workers	Number of threads for the extraction phase. Use 1 for deterministic sequential execution. TYPE: int
show_progress	When True, emit a tqdm progress bar over the region tasks. Requires tqdm to be installed; falls back to plain logging when absent. TYPE: bool

run

run(*, regions=None, variants=None)

Run the pipeline and return fully annotated variants.

Exactly one of regions or variants must be supplied.

• regions — each region is submitted as one extraction task.
• variants — each variant is converted to a point region (pos, pos) and then de-duplicated by region to avoid redundant store reads. Variants without a valid position are skipped.

PARAMETER	DESCRIPTION
regions	Genomic windows to extract variants from. TYPE: Sequence[GenomicPosition] | None DEFAULT: None
variants	Specific variants whose positions define point regions. TYPE: Sequence[Variant] | None DEFAULT: None

RETURNS	DESCRIPTION
list[VariantAnnotation]	Fully annotated list of
list[VariantAnnotation]	class:~iris.domain.genomics.entities.VariantAnnotation.

RAISES	DESCRIPTION
ValueError	If both or neither of regions and variants are supplied.

Source code in src/iris/adapters/driven/pipeline.py

def run(
    self,
    *,
    regions: Sequence[GenomicPosition] | None = None,
    variants: Sequence[Variant] | None = None,
) -> list[VariantAnnotation]:
    """Run the pipeline and return fully annotated variants.

    Exactly one of *regions* or *variants* must be supplied.

    * **regions** — each region is submitted as one extraction task.
    * **variants** — each variant is converted to a point region
      ``(pos, pos)`` and then de-duplicated by region to avoid redundant
      store reads.  Variants without a valid position are skipped.

    Args:
        regions: Genomic windows to extract variants from.
        variants: Specific variants whose positions define point regions.

    Returns:
        Fully annotated list of
        :class:`~iris.domain.genomics.entities.VariantAnnotation`.

    Raises:
        ValueError: If both or neither of *regions* and *variants* are
            supplied.
    """
    if (regions is None) == (variants is None):
        raise ValueError("Exactly one of 'regions' or 'variants' must be provided.")

    tasks = self._resolve_tasks(regions=regions, variants=variants)
    logger.info(
        "pipeline start — %s tasks · workers=%d · %s",
        f"{len(tasks):,}",
        self.workers,
        _fmt_rss(_rss_gb()),
    )
    extracted = self._extract_parallel(tasks)
    logger.info(
        "extraction done — %s variants · %s", f"{len(extracted):,}", _fmt_rss(_rss_gb())
    )
    return self._annotate_sequential(extracted)