RecGenes — Zarr WES + ClinVar + LoF¶

Variante del pipeline original donde el paso 2 cambia: las variantes se obtienen del Zarr WES (genotipos de todas las muestras) en lugar de los archivos tabix MCPS-AF.

Pasos:

1. BED → GenomicPositions — regiones génicas por cromosoma.
2. Zarr WES → VariantAnnotation + GenotypeCount — variantes dentro de las regiones; se calculan n_hom, n_het, n_hom_ref, n_missing para todas las muestras usando operaciones numpy vectorizadas.
3. ClinVar — significancia clínica.
4. WES LoF — transcritos MANE Select + score LoF.
5. Estadísticas de población — número de homocigotos por variante y mutaciones promedio por participante.

In [ ]:

%load_ext autoreload
%autoreload 2

%load_ext autoreload %autoreload 2

In [ ]:

from pathlib import Path

BED_DIR      = Path("/mnt/cephfs/orgs/home/gen/repos/recgenes/genes_regions")
ZARR_DIR     = Path("/mnt/cephfs/hot_nvme/mcps/whole_exome_sequencing/oqfe-svm/zar_files")
CLINVAR_TSV  = Path("/mnt/cephfs/hot_nvme/clinvar/GRCh38/clinvar.tsv.gz")
WES_LOF_TABIX = Path(
    "/mnt/cephfs/hot_nvme/mcps/whole_exome_sequencing/oqfe-svm"
    "/wes-svm-filtered.sites.norm.vep.loftee.tabix.tsv.gz"
)
WORKERS = 4

from pathlib import Path BED_DIR = Path("/mnt/cephfs/orgs/home/gen/repos/recgenes/genes_regions") ZARR_DIR = Path("/mnt/cephfs/hot_nvme/mcps/whole_exome_sequencing/oqfe-svm/zar_files") CLINVAR_TSV = Path("/mnt/cephfs/hot_nvme/clinvar/GRCh38/clinvar.tsv.gz") WES_LOF_TABIX = Path( "/mnt/cephfs/hot_nvme/mcps/whole_exome_sequencing/oqfe-svm" "/wes-svm-filtered.sites.norm.vep.loftee.tabix.tsv.gz" ) WORKERS = 4

In [ ]:

from __future__ import annotations

from bisect import bisect_right
from pathlib import Path

import numpy as np
import zarr
import attrs
from attrs import define

from iris.adapters.driven.annotators.clinvar import ClinvarAnnotator
from iris.adapters.driven.annotators.tabix import TabixVariantAnnotator
from iris.adapters.driven.readers.bed import BedRegionReader
from iris.adapters.driven.readers.tabix import _Key
from iris.domain.generic.abstracts import Metadata
from iris.domain.genomics.entities import (
    Gene,
    Transcript,
    TranscriptVariant,
    Variant,
    VariantAnnotation,
)
from iris.domain.genomics.objects import (
    Allele,
    Biotype,
    GenotypeCount,
    GenomicPosition,
    ScorePredictionRank,
    TranscriptNomenclature,
    VariantNomenclature,
    VariantScore,
    VariantScoreType,
    VariantType,
)
from iris.domain.genomics.utils import normalize_chromosome

print("Imports OK")

from __future__ import annotations from bisect import bisect_right from pathlib import Path import numpy as np import zarr import attrs from attrs import define from iris.adapters.driven.annotators.clinvar import ClinvarAnnotator from iris.adapters.driven.annotators.tabix import TabixVariantAnnotator from iris.adapters.driven.readers.bed import BedRegionReader from iris.adapters.driven.readers.tabix import _Key from iris.domain.generic.abstracts import Metadata from iris.domain.genomics.entities import ( Gene, Transcript, TranscriptVariant, Variant, VariantAnnotation, ) from iris.domain.genomics.objects import ( Allele, Biotype, GenotypeCount, GenomicPosition, ScorePredictionRank, TranscriptNomenclature, VariantNomenclature, VariantScore, VariantScoreType, VariantType, ) from iris.domain.genomics.utils import normalize_chromosome print("Imports OK")

---

1. Leer regiones BED por cromosoma¶

In [ ]:

bed_reader = BedRegionReader()
regions_by_chrom: dict[str, list[GenomicPosition]] = {}

for bed_file in sorted(BED_DIR.glob("chr*.bed")):
    with open(bed_file) as fh:
        regions = bed_reader.read(fh, assembly="GRCh38")
    if regions:
        chrom = regions[0].chromosome
        regions_by_chrom[chrom] = regions

total_regions = sum(len(v) for v in regions_by_chrom.values())
print(f"{len(regions_by_chrom)} cromosomas · {total_regions} regiones totales")

bed_reader = BedRegionReader() regions_by_chrom: dict[str, list[GenomicPosition]] = {} for bed_file in sorted(BED_DIR.glob("chr*.bed")): with open(bed_file) as fh: regions = bed_reader.read(fh, assembly="GRCh38") if regions: chrom = regions[0].chromosome regions_by_chrom[chrom] = regions total_regions = sum(len(v) for v in regions_by_chrom.values()) print(f"{len(regions_by_chrom)} cromosomas · {total_regions} regiones totales")

In [ ]:

_RegionIndex = list[tuple[int, int, str]]


def _build_gene_index(regions: list[GenomicPosition]) -> _RegionIndex:
    return sorted(
        [
            (r.start, r.end, r.name)
            for r in regions
            if r.start is not None and r.end is not None and r.name
        ],
        key=lambda t: t[0],
    )


def _find_gene(pos: int, index: _RegionIndex) -> str | None:
    starts = [t[0] for t in index]
    i = bisect_right(starts, pos) - 1
    if i < 0:
        return None
    start, end, name = index[i]
    return name if start <= pos <= end else None


gene_index_by_chrom: dict[str, _RegionIndex] = {
    chrom: _build_gene_index(regs) for chrom, regs in regions_by_chrom.items()
}
print("Índice de genes construido.")

_RegionIndex = list[tuple[int, int, str]] def _build_gene_index(regions: list[GenomicPosition]) -> _RegionIndex: return sorted( [ (r.start, r.end, r.name) for r in regions if r.start is not None and r.end is not None and r.name ], key=lambda t: t[0], ) def _find_gene(pos: int, index: _RegionIndex) -> str | None: starts = [t[0] for t in index] i = bisect_right(starts, pos) - 1 if i < 0: return None start, end, name = index[i] return name if start <= pos <= end else None gene_index_by_chrom: dict[str, _RegionIndex] = { chrom: _build_gene_index(regs) for chrom, regs in regions_by_chrom.items() } print("Índice de genes construido.")

---

2. Zarr WES → VariantAnnotation + GenotypeCount¶

Reemplaza el paso 2 original (tabix MCPS-AF).
Para cada cromosoma:

1. Abrir el Zarr del cromosoma.
2. Construir una máscara booleana sobre (variant_contig, variant_position) para las regiones BED — todo con numpy, sin iterar fila a fila.
3. Cargar call_genotype[mask, :, :] — solo las variantes en región, todos los samples.
4. Calcular GenotypeCount vectorizado: n_hom_ref, n_het, n_hom_alt, n_missing para cada variante filtrada en una sola pasada.
5. Construir VariantAnnotation con AlleleFrequency derivada de GenotypeCount y los counts brutos en Metadata.

In [ ]:

def _open_zarr(zarr_dir: Path, chrom: str):
    """Abre el Zarr del cromosoma; prueba varios patrones de nombre."""
    for name in [f"chr{chrom}.zarr", f"chr{chrom}", chrom]:
        p = zarr_dir / name
        if p.exists():
            return zarr.open(str(p), mode="r")
    return None


def _region_variant_mask(
    v_contig: np.ndarray,
    v_position: np.ndarray,
    chrom_idx: int,
    regions: list[GenomicPosition],
) -> np.ndarray:
    """Máscara booleana (n_variants,) para variantes en las regiones del cromosoma.

    v_position es 0-based (vcf2zarr default); las regiones BED están en
    1-based inclusive → se convierten aquí.
    """
    on_chrom = v_contig == chrom_idx
    in_region = np.zeros(len(v_contig), dtype=bool)
    for r in regions:
        if r.start is None or r.end is None:
            continue
        start_0 = r.start - 1          # 1-based → 0-based
        end_0   = r.end   - 1          # inclusive
        in_region |= on_chrom & (v_position >= start_0) & (v_position <= end_0)
    return in_region


def _genotype_counts_vectorized(
    gt: np.ndarray,           # (n_variants, n_samples, ploidy)
    gt_mask: np.ndarray | None,  # (n_variants, n_samples, ploidy) bool
    alt_idx: int = 1,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Calcula n_hom_ref, n_het, n_hom_alt, n_missing por variante (eje 0).

    Retorna cuatro arrays de forma (n_variants,).
    """
    if gt_mask is not None:
        missing_allele = gt_mask.any(axis=-1)     # (n_variants, n_samples)
        gt = np.where(gt_mask, -1, gt)
    else:
        missing_allele = (gt < 0).any(axis=-1)

    called = ~missing_allele                       # (n_variants, n_samples)
    has_alt  = (gt == alt_idx).any(axis=-1)        # (n_variants, n_samples)
    all_alt  = (gt == alt_idx).all(axis=-1)        # (n_variants, n_samples)
    all_ref  = (gt == 0).all(axis=-1)              # (n_variants, n_samples)

    n_hom_alt = (called & all_alt).sum(axis=1).astype(int)
    n_het     = (called & has_alt & ~all_alt).sum(axis=1).astype(int)
    n_hom_ref = (called & all_ref).sum(axis=1).astype(int)
    n_missing = (~called).sum(axis=1).astype(int)

    return n_hom_ref, n_het, n_hom_alt, n_missing


print("Helpers definidos.")

def _open_zarr(zarr_dir: Path, chrom: str): """Abre el Zarr del cromosoma; prueba varios patrones de nombre.""" for name in [f"chr{chrom}.zarr", f"chr{chrom}", chrom]: p = zarr_dir / name if p.exists(): return zarr.open(str(p), mode="r") return None def _region_variant_mask( v_contig: np.ndarray, v_position: np.ndarray, chrom_idx: int, regions: list[GenomicPosition], ) -> np.ndarray: """Máscara booleana (n_variants,) para variantes en las regiones del cromosoma. v_position es 0-based (vcf2zarr default); las regiones BED están en 1-based inclusive → se convierten aquí. """ on_chrom = v_contig == chrom_idx in_region = np.zeros(len(v_contig), dtype=bool) for r in regions: if r.start is None or r.end is None: continue start_0 = r.start - 1 # 1-based → 0-based end_0 = r.end - 1 # inclusive in_region |= on_chrom & (v_position >= start_0) & (v_position <= end_0) return in_region def _genotype_counts_vectorized( gt: np.ndarray, # (n_variants, n_samples, ploidy) gt_mask: np.ndarray | None, # (n_variants, n_samples, ploidy) bool alt_idx: int = 1, ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """Calcula n_hom_ref, n_het, n_hom_alt, n_missing por variante (eje 0). Retorna cuatro arrays de forma (n_variants,). """ if gt_mask is not None: missing_allele = gt_mask.any(axis=-1) # (n_variants, n_samples) gt = np.where(gt_mask, -1, gt) else: missing_allele = (gt < 0).any(axis=-1) called = ~missing_allele # (n_variants, n_samples) has_alt = (gt == alt_idx).any(axis=-1) # (n_variants, n_samples) all_alt = (gt == alt_idx).all(axis=-1) # (n_variants, n_samples) all_ref = (gt == 0).all(axis=-1) # (n_variants, n_samples) n_hom_alt = (called & all_alt).sum(axis=1).astype(int) n_het = (called & has_alt & ~all_alt).sum(axis=1).astype(int) n_hom_ref = (called & all_ref).sum(axis=1).astype(int) n_missing = (~called).sum(axis=1).astype(int) return n_hom_ref, n_het, n_hom_alt, n_missing print("Helpers definidos.")

In [ ]:

all_variants: list[VariantAnnotation] = []
total_samples: int | None = None

for chrom, chrom_regions in sorted(regions_by_chrom.items()):
    store = _open_zarr(ZARR_DIR, chrom)
    if store is None:
        print(f"  chr{chrom}: Zarr no encontrado, saltando.")
        continue

    # ── Lookup arrays (pequeños, caben en RAM) ──────────────────────────
    contig_id  = store["contig_id"][:]
    sample_ids = list(store["sample_id"][:]) if "sample_id" in store else []
    n_samples  = len(sample_ids)
    if total_samples is None:
        total_samples = n_samples

    # Índice del cromosoma en el Zarr
    chrom_norm = normalize_chromosome(f"chr{chrom}")
    chrom_indices = [
        i for i, c in enumerate(contig_id)
        if normalize_chromosome(str(c)) == chrom_norm
    ]
    if not chrom_indices:
        print(f"  chr{chrom}: cromosoma no encontrado en el Zarr.")
        continue
    chrom_idx = chrom_indices[0]

    # ── Máscara de región ───────────────────────────────────────────────
    v_contig   = store["variant_contig"][:]
    v_position = store["variant_position"][:]
    mask       = _region_variant_mask(v_contig, v_position, chrom_idx, chrom_regions)
    indices    = np.where(mask)[0]

    if len(indices) == 0:
        print(f"  chr{chrom}: 0 variantes en región.")
        continue

    # ── Arrays de variante para las posiciones filtradas ────────────────
    v_allele = store["variant_allele"][indices]
    v_id     = store["variant_id"][indices] if "variant_id" in store else None

    # ── Genotipos: cargar solo las filas filtradas ───────────────────────
    # call_genotype shape: (n_variants_total, n_samples, ploidy)
    gt_full  = store["call_genotype"][indices]           # (n_filt, n_samples, ploidy)
    gt_mask_ = None
    if "call_genotype_mask" in store:
        gt_mask_ = store["call_genotype_mask"][indices]   # (n_filt, n_samples, ploidy)

    gene_index = gene_index_by_chrom[chrom]
    chrom_variants: list[VariantAnnotation] = []

    for local_i, global_i in enumerate(indices):
        pos_0  = int(v_position[global_i])
        pos_1  = pos_0 + 1                               # 0-based → 1-based
        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

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

        # Genotipo de esta variante: (n_samples, ploidy)
        gt_row   = gt_full[local_i]                      # (n_samples, ploidy)
        mask_row = gt_mask_[local_i] if gt_mask_ is not None else None

        for alt_idx, alt in enumerate(alts, start=1):
            # Calcular GenotypeCount para este alelo
            gt_2d   = gt_row[np.newaxis]                  # (1, n_samples, ploidy)
            mask_2d = mask_row[np.newaxis] if mask_row is not None else None
            n_hom_ref, n_het, n_hom_alt, n_missing = _genotype_counts_vectorized(
                gt_2d, mask_2d, alt_idx=alt_idx
            )
            gc = GenotypeCount(
                n_hom_ref=int(n_hom_ref[0]),
                n_het=int(n_het[0]),
                n_hom_alt=int(n_hom_alt[0]),
                n_missing=int(n_missing[0]),
                allele=alt,
            )

            key = (chrom_norm, pos_1, ref, alt)
            gene_name = _find_gene(pos_1, gene_index)
            genes = (Gene.from_symbol(gene_name),) if gene_name else ()

            chrom_variants.append(
                VariantAnnotation(
                    variant=Variant(
                        variant_type=VariantType.from_key(key),
                        nomenclature=VariantNomenclature.from_key(
                            key,
                            dbsnp_id=vid_raw if vid_raw and vid_raw.startswith("rs") else None,
                        ),
                        alleles=Allele(reference=ref, alternate=(alt,)),
                        position=GenomicPosition(
                            chromosome=chrom_norm,
                            assembly="GRCh38",
                            start=pos_1,
                            end=pos_1 + len(ref) - 1,
                        ),
                    ),
                    genes=genes,
                    allele_frequencies=(gc.to_allele_frequency("zarr_wes"),),
                    metadata=Metadata(
                        source="zarr_wes",
                        data={
                            "n_hom":    gc.n_hom_alt,
                            "n_het":    gc.n_het,
                            "n_hom_ref":gc.n_hom_ref,
                            "n_missing":gc.n_missing,
                            "n_samples":n_samples,
                            "ac":       gc.ac,
                            "an":       gc.an,
                        },
                    ),
                )
            )

    all_variants.extend(chrom_variants)
    print(f"  chr{chrom}: {len(chrom_variants)} variantes ({n_samples} muestras)")

print(f"\nTotal: {len(all_variants)} variantes · {total_samples} muestras")

all_variants: list[VariantAnnotation] = [] total_samples: int | None = None for chrom, chrom_regions in sorted(regions_by_chrom.items()): store = _open_zarr(ZARR_DIR, chrom) if store is None: print(f" chr{chrom}: Zarr no encontrado, saltando.") continue # ── Lookup arrays (pequeños, caben en RAM) ────────────────────────── contig_id = store["contig_id"][:] sample_ids = list(store["sample_id"][:]) if "sample_id" in store else [] n_samples = len(sample_ids) if total_samples is None: total_samples = n_samples # Índice del cromosoma en el Zarr chrom_norm = normalize_chromosome(f"chr{chrom}") chrom_indices = [ i for i, c in enumerate(contig_id) if normalize_chromosome(str(c)) == chrom_norm ] if not chrom_indices: print(f" chr{chrom}: cromosoma no encontrado en el Zarr.") continue chrom_idx = chrom_indices[0] # ── Máscara de región ─────────────────────────────────────────────── v_contig = store["variant_contig"][:] v_position = store["variant_position"][:] mask = _region_variant_mask(v_contig, v_position, chrom_idx, chrom_regions) indices = np.where(mask)[0] if len(indices) == 0: print(f" chr{chrom}: 0 variantes en región.") continue # ── Arrays de variante para las posiciones filtradas ──────────────── v_allele = store["variant_allele"][indices] v_id = store["variant_id"][indices] if "variant_id" in store else None # ── Genotipos: cargar solo las filas filtradas ─────────────────────── # call_genotype shape: (n_variants_total, n_samples, ploidy) gt_full = store["call_genotype"][indices] # (n_filt, n_samples, ploidy) gt_mask_ = None if "call_genotype_mask" in store: gt_mask_ = store["call_genotype_mask"][indices] # (n_filt, n_samples, ploidy) gene_index = gene_index_by_chrom[chrom] chrom_variants: list[VariantAnnotation] = [] for local_i, global_i in enumerate(indices): pos_0 = int(v_position[global_i]) pos_1 = pos_0 + 1 # 0-based → 1-based 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 vid_raw = str(v_id[local_i]) if v_id is not None else None if vid_raw in ("", "."): vid_raw = None # Genotipo de esta variante: (n_samples, ploidy) gt_row = gt_full[local_i] # (n_samples, ploidy) mask_row = gt_mask_[local_i] if gt_mask_ is not None else None for alt_idx, alt in enumerate(alts, start=1): # Calcular GenotypeCount para este alelo gt_2d = gt_row[np.newaxis] # (1, n_samples, ploidy) mask_2d = mask_row[np.newaxis] if mask_row is not None else None n_hom_ref, n_het, n_hom_alt, n_missing = _genotype_counts_vectorized( gt_2d, mask_2d, alt_idx=alt_idx ) gc = GenotypeCount( n_hom_ref=int(n_hom_ref[0]), n_het=int(n_het[0]), n_hom_alt=int(n_hom_alt[0]), n_missing=int(n_missing[0]), allele=alt, ) key = (chrom_norm, pos_1, ref, alt) gene_name = _find_gene(pos_1, gene_index) genes = (Gene.from_symbol(gene_name),) if gene_name else () chrom_variants.append( VariantAnnotation( variant=Variant( variant_type=VariantType.from_key(key), nomenclature=VariantNomenclature.from_key( key, dbsnp_id=vid_raw if vid_raw and vid_raw.startswith("rs") else None, ), alleles=Allele(reference=ref, alternate=(alt,)), position=GenomicPosition( chromosome=chrom_norm, assembly="GRCh38", start=pos_1, end=pos_1 + len(ref) - 1, ), ), genes=genes, allele_frequencies=(gc.to_allele_frequency("zarr_wes"),), metadata=Metadata( source="zarr_wes", data={ "n_hom": gc.n_hom_alt, "n_het": gc.n_het, "n_hom_ref":gc.n_hom_ref, "n_missing":gc.n_missing, "n_samples":n_samples, "ac": gc.ac, "an": gc.an, }, ), ) ) all_variants.extend(chrom_variants) print(f" chr{chrom}: {len(chrom_variants)} variantes ({n_samples} muestras)") print(f"\nTotal: {len(all_variants)} variantes · {total_samples} muestras")

---

3. Anotar significancia clínica (ClinVar)¶

In [ ]:

clinvar = ClinvarAnnotator(path=CLINVAR_TSV)
all_variants = clinvar.annotate(all_variants)
print(f"ClinVar anotado: {len(all_variants)} variantes")

clinvar = ClinvarAnnotator(path=CLINVAR_TSV) all_variants = clinvar.annotate(all_variants) print(f"ClinVar anotado: {len(all_variants)} variantes")

In [ ]:

def _keep_most_pathogenic(va: VariantAnnotation) -> VariantAnnotation:
    if len(va.classifications) <= 1:
        return va
    best = max(va.classifications, key=lambda c: c.classification.pathogenicity_rank)
    return attrs.evolve(va, classifications=(best,))


all_variants = [_keep_most_pathogenic(va) for va in all_variants]
n_classified = sum(1 for va in all_variants if va.classifications)
print(f"Variantes con clasificación ClinVar: {n_classified} / {len(all_variants)}")

def _keep_most_pathogenic(va: VariantAnnotation) -> VariantAnnotation: if len(va.classifications) <= 1: return va best = max(va.classifications, key=lambda c: c.classification.pathogenicity_rank) return attrs.evolve(va, classifications=(best,)) all_variants = [_keep_most_pathogenic(va) for va in all_variants] n_classified = sum(1 for va in all_variants if va.classifications) print(f"Variantes con clasificación ClinVar: {n_classified} / {len(all_variants)}")

---

4. Anotar transcritos + LoF (WES VEP/LOFTEE)¶

In [ ]:

@define
class WesLofteeAnnotator(TabixVariantAnnotator):
    path: str | Path = ""

    def annotate(
        self,
        items: list[VariantAnnotation],
        *,
        workers: int = 1,
    ) -> list[VariantAnnotation]:
        return self._annotate_rows(
            items,
            Path(self.path),
            key_fn=self._key_fn,
            enrich_fn=self._enrich_fn,
            columns=["Uploaded_variation", "Feature", "MANE_SELECT", "Consequence", "LoF"],
            workers=workers,
        )

    @staticmethod
    def _key_fn(row: dict[str, str]) -> _Key | None:
        parts = row.get("Uploaded_variation", "").split(":")
        if len(parts) != 4:
            return None
        chrom, pos_s, ref, alt = parts
        try:
            return (normalize_chromosome(chrom), int(pos_s), ref, alt)
        except ValueError:
            return None

    @staticmethod
    def _enrich_fn(va: VariantAnnotation, row: dict[str, str]) -> VariantAnnotation:
        feature  = row.get("Feature", "-")
        mane_id  = row.get("MANE_SELECT", "-")
        consq    = row.get("Consequence", "")
        lof      = row.get("LoF", "-")
        updated  = va
        tx_name  = mane_id if (mane_id and mane_id != "-") else feature
        is_mane  = bool(mane_id and mane_id != "-")
        if consq and consq != "-" and tx_name and tx_name != "-":
            tx = Transcript(
                biotype=Biotype.PROTEIN_CODING,
                nomenclature=TranscriptNomenclature.from_name(tx_name),
                is_mane_select=True if is_mane else None,
            )
            tc = TranscriptVariant(
                transcript=tx, consequences=tuple(consq.split("&"))
            )
            updated = attrs.evolve(updated, transcripts=(*updated.transcripts, tc))
        if lof and lof != "-" and not any(s.name == "LoF" for s in updated.scores):
            updated = attrs.evolve(
                updated,
                scores=(
                    *updated.scores,
                    VariantScore(
                        name="LoF",
                        type=VariantScoreType.PATHOGENICITY,
                        value=None,
                        prediction=ScorePredictionRank.VERY_HIGH,
                    ),
                ),
            )
        return updated

@define class WesLofteeAnnotator(TabixVariantAnnotator): path: str | Path = "" def annotate( self, items: list[VariantAnnotation], *, workers: int = 1, ) -> list[VariantAnnotation]: return self._annotate_rows( items, Path(self.path), key_fn=self._key_fn, enrich_fn=self._enrich_fn, columns=["Uploaded_variation", "Feature", "MANE_SELECT", "Consequence", "LoF"], workers=workers, ) @staticmethod def _key_fn(row: dict[str, str]) -> _Key | None: parts = row.get("Uploaded_variation", "").split(":") if len(parts) != 4: return None chrom, pos_s, ref, alt = parts try: return (normalize_chromosome(chrom), int(pos_s), ref, alt) except ValueError: return None @staticmethod def _enrich_fn(va: VariantAnnotation, row: dict[str, str]) -> VariantAnnotation: feature = row.get("Feature", "-") mane_id = row.get("MANE_SELECT", "-") consq = row.get("Consequence", "") lof = row.get("LoF", "-") updated = va tx_name = mane_id if (mane_id and mane_id != "-") else feature is_mane = bool(mane_id and mane_id != "-") if consq and consq != "-" and tx_name and tx_name != "-": tx = Transcript( biotype=Biotype.PROTEIN_CODING, nomenclature=TranscriptNomenclature.from_name(tx_name), is_mane_select=True if is_mane else None, ) tc = TranscriptVariant( transcript=tx, consequences=tuple(consq.split("&")) ) updated = attrs.evolve(updated, transcripts=(*updated.transcripts, tc)) if lof and lof != "-" and not any(s.name == "LoF" for s in updated.scores): updated = attrs.evolve( updated, scores=( *updated.scores, VariantScore( name="LoF", type=VariantScoreType.PATHOGENICITY, value=None, prediction=ScorePredictionRank.VERY_HIGH, ), ), ) return updated

In [ ]:

loftee = WesLofteeAnnotator(path=WES_LOF_TABIX)
all_variants = loftee.annotate(all_variants, workers=WORKERS)

n_lof = sum(1 for va in all_variants if any(s.name == "LoF" for s in va.scores))
print(f"Variantes con score LoF: {n_lof} / {len(all_variants)}")

loftee = WesLofteeAnnotator(path=WES_LOF_TABIX) all_variants = loftee.annotate(all_variants, workers=WORKERS) n_lof = sum(1 for va in all_variants if any(s.name == "LoF" for s in va.scores)) print(f"Variantes con score LoF: {n_lof} / {len(all_variants)}")

---

5. Estadísticas de población¶

5.1 Distribución de homocigotos por variante¶

Cuántas muestras tienen el genotipo homocigoto para cada variante (n_hom).

5.2 Mutaciones promedio por participante¶

Para cada muestra: cuántas variantes de la lista lleva (≥1 alelo ALT).
El promedio sobre las N muestras es:

$$\bar{m} = \frac{\sum_v (n_{\text{het},v} + n_{\text{hom},v})}{N_{\text{muestras}}}$$

In [ ]:

from collections import Counter

# ── Distribución de n_hom ───────────────────────────────────────────────────
n_hom_values = [
    va.metadata.get("n_hom", 0)
    for va in all_variants
    if va.metadata is not None
]

hom_dist = Counter(n_hom_values)
n_with_hom = sum(1 for v in n_hom_values if v > 0)
max_hom    = max(n_hom_values) if n_hom_values else 0

print(f"Variantes con ≥1 homocigoto : {n_with_hom:,}")
print(f"Máx. homocigotos en una variante: {max_hom:,}")
print()

# Top 20 variantes por n_hom
print(f"{'Variante':<35} {'Gen':<12} {'n_hom':>8} {'n_het':>8} {'AF':>8}")
print("-" * 73)
top_hom = sorted(
    all_variants,
    key=lambda va: va.metadata.get("n_hom", 0) if va.metadata else 0,
    reverse=True,
)[:20]

for va in top_hom:
    n_hom  = va.metadata.get("n_hom", 0) if va.metadata else 0
    n_het  = va.metadata.get("n_het", 0) if va.metadata else 0
    gene   = va.genes[0].symbol if va.genes else "-"
    af_obj = va.af_for("zarr_wes")
    af_str = f"{af_obj.frequency:.4f}" if af_obj and af_obj.frequency else "-"
    print(f"{va.variant.cpra:<35} {gene:<12} {n_hom:>8,} {n_het:>8,} {af_str:>8}")

from collections import Counter # ── Distribución de n_hom ─────────────────────────────────────────────────── n_hom_values = [ va.metadata.get("n_hom", 0) for va in all_variants if va.metadata is not None ] hom_dist = Counter(n_hom_values) n_with_hom = sum(1 for v in n_hom_values if v > 0) max_hom = max(n_hom_values) if n_hom_values else 0 print(f"Variantes con ≥1 homocigoto : {n_with_hom:,}") print(f"Máx. homocigotos en una variante: {max_hom:,}") print() # Top 20 variantes por n_hom print(f"{'Variante':<35} {'Gen':<12} {'n_hom':>8} {'n_het':>8} {'AF':>8}") print("-" * 73) top_hom = sorted( all_variants, key=lambda va: va.metadata.get("n_hom", 0) if va.metadata else 0, reverse=True, )[:20] for va in top_hom: n_hom = va.metadata.get("n_hom", 0) if va.metadata else 0 n_het = va.metadata.get("n_het", 0) if va.metadata else 0 gene = va.genes[0].symbol if va.genes else "-" af_obj = va.af_for("zarr_wes") af_str = f"{af_obj.frequency:.4f}" if af_obj and af_obj.frequency else "-" print(f"{va.variant.cpra:<35} {gene:<12} {n_hom:>8,} {n_het:>8,} {af_str:>8}")

In [ ]:

# ── Mutaciones promedio por participante ────────────────────────────────────
#
# total_carriers = Σ_v (n_het_v + n_hom_v)  (portadores de ≥1 ALT en la variante v)
# avg_mutations  = total_carriers / N_muestras
#
# Interpretación: número esperado de variantes de este panel que
# lleva un participante aleatorio del WES.

if total_samples and total_samples > 0:
    total_carriers = sum(
        (va.metadata.get("n_het", 0) + va.metadata.get("n_hom", 0))
        for va in all_variants
        if va.metadata is not None
    )
    avg_mutations = total_carriers / total_samples

    print(f"Muestras totales          : {total_samples:,}")
    print(f"Total portadores (v×s)    : {total_carriers:,}")
    print(f"Mutaciones promedio/part. : {avg_mutations:.4f}")

    # Desglose por gen
    from collections import defaultdict
    carriers_by_gene: dict[str, int] = defaultdict(int)
    for va in all_variants:
        if va.metadata is None:
            continue
        gene = va.genes[0].symbol if va.genes else "DESCONOCIDO"
        carriers_by_gene[gene] += (
            va.metadata.get("n_het", 0) + va.metadata.get("n_hom", 0)
        )

    print()
    print(f"{'Gen':<15} {'Portadores':>12} {'Avg/part':>10}")
    print("-" * 40)
    for gene, cnt in sorted(carriers_by_gene.items(), key=lambda x: -x[1])[:20]:
        print(f"{gene:<15} {cnt:>12,} {cnt/total_samples:>10.4f}")

# ── Mutaciones promedio por participante ──────────────────────────────────── # # total_carriers = Σ_v (n_het_v + n_hom_v) (portadores de ≥1 ALT en la variante v) # avg_mutations = total_carriers / N_muestras # # Interpretación: número esperado de variantes de este panel que # lleva un participante aleatorio del WES. if total_samples and total_samples > 0: total_carriers = sum( (va.metadata.get("n_het", 0) + va.metadata.get("n_hom", 0)) for va in all_variants if va.metadata is not None ) avg_mutations = total_carriers / total_samples print(f"Muestras totales : {total_samples:,}") print(f"Total portadores (v×s) : {total_carriers:,}") print(f"Mutaciones promedio/part. : {avg_mutations:.4f}") # Desglose por gen from collections import defaultdict carriers_by_gene: dict[str, int] = defaultdict(int) for va in all_variants: if va.metadata is None: continue gene = va.genes[0].symbol if va.genes else "DESCONOCIDO" carriers_by_gene[gene] += ( va.metadata.get("n_het", 0) + va.metadata.get("n_hom", 0) ) print() print(f"{'Gen':<15} {'Portadores':>12} {'Avg/part':>10}") print("-" * 40) for gene, cnt in sorted(carriers_by_gene.items(), key=lambda x: -x[1])[:20]: print(f"{gene:<15} {cnt:>12,} {cnt/total_samples:>10.4f}")

---

6. Exportar TSV¶

In [ ]:

import csv

OUTPUT_PATH = Path("recgenes_zarr_variants.tsv")

_HEADERS = [
    "variant_id", "chromosome", "start", "end", "ref", "alt",
    "variant_type", "genes", "consequence",
    "clinvar_sig", "clinvar_evidence",
    "zarr_af", "zarr_ac", "zarr_an",
    "n_hom", "n_het", "n_hom_ref", "n_missing", "n_samples",
    "lof",
]


def _get_af(va: VariantAnnotation, population: str) -> str:
    af = va.af_for(population)
    return str(af.frequency) if af and af.frequency is not None else ""


def _meta(va: VariantAnnotation, key: str) -> str:
    if va.metadata is None:
        return ""
    v = va.metadata.get(key)
    return str(v) if v is not None else ""


def _export_row(va: VariantAnnotation) -> list:
    v    = va.variant
    pos  = v.position
    al   = v.alleles
    clas = va.top_classification
    lof  = next((s.prediction.name for s in va.scores if s.name == "LoF" and s.prediction), "")
    af   = va.af_for("zarr_wes")
    return [
        v.cpra,
        pos.chromosome if pos else "",
        pos.start if pos else "",
        pos.end if pos else "",
        al.reference if al else "",
        al.alt if al else "",
        v.variant_type.name,
        ";".join(va.gene_symbols),
        va.most_severe_consequence or "",
        clas.classification.name if clas else "",
        clas.evidence.name if clas and clas.evidence else "",
        str(af.frequency) if af and af.frequency is not None else "",
        str(af.count) if af and af.count is not None else "",
        str(af.total) if af and af.total is not None else "",
        _meta(va, "n_hom"),
        _meta(va, "n_het"),
        _meta(va, "n_hom_ref"),
        _meta(va, "n_missing"),
        _meta(va, "n_samples"),
        lof,
    ]


with OUTPUT_PATH.open("w", newline="") as fh:
    writer = csv.writer(fh, delimiter="\t")
    writer.writerow(_HEADERS)
    for va in all_variants:
        writer.writerow(_export_row(va))

print(f"Exportadas {len(all_variants):,} variantes → {OUTPUT_PATH}")

import csv OUTPUT_PATH = Path("recgenes_zarr_variants.tsv") _HEADERS = [ "variant_id", "chromosome", "start", "end", "ref", "alt", "variant_type", "genes", "consequence", "clinvar_sig", "clinvar_evidence", "zarr_af", "zarr_ac", "zarr_an", "n_hom", "n_het", "n_hom_ref", "n_missing", "n_samples", "lof", ] def _get_af(va: VariantAnnotation, population: str) -> str: af = va.af_for(population) return str(af.frequency) if af and af.frequency is not None else "" def _meta(va: VariantAnnotation, key: str) -> str: if va.metadata is None: return "" v = va.metadata.get(key) return str(v) if v is not None else "" def _export_row(va: VariantAnnotation) -> list: v = va.variant pos = v.position al = v.alleles clas = va.top_classification lof = next((s.prediction.name for s in va.scores if s.name == "LoF" and s.prediction), "") af = va.af_for("zarr_wes") return [ v.cpra, pos.chromosome if pos else "", pos.start if pos else "", pos.end if pos else "", al.reference if al else "", al.alt if al else "", v.variant_type.name, ";".join(va.gene_symbols), va.most_severe_consequence or "", clas.classification.name if clas else "", clas.evidence.name if clas and clas.evidence else "", str(af.frequency) if af and af.frequency is not None else "", str(af.count) if af and af.count is not None else "", str(af.total) if af and af.total is not None else "", _meta(va, "n_hom"), _meta(va, "n_het"), _meta(va, "n_hom_ref"), _meta(va, "n_missing"), _meta(va, "n_samples"), lof, ] with OUTPUT_PATH.open("w", newline="") as fh: writer = csv.writer(fh, delimiter="\t") writer.writerow(_HEADERS) for va in all_variants: writer.writerow(_export_row(va)) print(f"Exportadas {len(all_variants):,} variantes → {OUTPUT_PATH}")