diff --git a/Makefile b/Makefile index 9dfcfd5..34c74a2 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -.PHONY: up down clean migrate test lint data loader pipeline annotate images kind serverless-deploy serverless-destroy gcp-configure gcp-secrets +.PHONY: up down clean migrate test lint data hpo demo-case loader pipeline annotate images kind serverless-deploy serverless-destroy gcp-configure gcp-secrets VCF ?= data/example.vcf.gz TAG ?= latest @@ -30,6 +30,12 @@ lint: data: ## download the public demo slice: GIAB HG002 + ClinVar, chr22 (see docs/data.md) scripts/fetch-demo-data.sh +hpo: ## load HPO gene-to-phenotype annotations, which the ranking matches against + scripts/load-hpo.py + +demo-case: ## build the simulated proband: GIAB background + one ClinVar pathogenic variant + scripts/make-demo-case.sh + loader: docker build -t rarelens/loader:dev -f pipeline/loader.Dockerfile pipeline diff --git a/README.md b/README.md index df02be9..f12913e 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,11 @@ # rarelens A small, end-to-end variant interpretation platform for rare genetic disease research. -Scientists upload a VCF, a Nextflow workflow annotates it with Ensembl VEP, a machine -learning model scores each variant, and results are browsable in a web app. +A case is a proband: a VCF plus the patient's phenotype (HPO terms). A Nextflow workflow +annotates the variants with Ensembl VEP, a model scores each one, and the app narrows +thousands of variants to a handful of candidates ranked against that phenotype — each +carrying the evidence for its rank, and each able to be shortlisted or dismissed with a +reason that ends up in a case report. This repository is a **self-training lab**. It exists so that one engineer can learn, in public, how a modern life-sciences platform is built end to end: full-stack application, @@ -28,11 +31,15 @@ clinical tool and makes no diagnostic claims. ```bash make up # postgres + api + web + mlflow via docker-compose make migrate # alembic upgrade head -make data # real public data: GIAB HG002 + ClinVar, chr22 (needs bcftools) +make hpo # HPO gene-to-phenotype annotations: what the ranking matches against +make demo-case # a simulated proband: GIAB background + one ClinVar pathogenic variant make test # api, ml, loader and web tests (no Docker needed for the DB tests) ``` -Then open http://localhost:5173. +Then open http://localhost:5173, create a case pointing at `data/proband-simulated.vcf.gz`, +give it the phenotype of the planted disease (for the default NF2 case: bilateral vestibular +schwannoma, sensorineural hearing impairment, tinnitus, meningioma, cataract), and analyse it. +The planted variant should come back ranked first. The docker-compose API has no Nextflow, so "Run VEP annotation" marks the job failed with the command to run instead. With Nextflow and Docker on the host, a VEP cache in `pipeline/cache/vep` diff --git a/api/alembic/versions/9a1c2d3e4f50_cases_phenotypes_decisions.py b/api/alembic/versions/9a1c2d3e4f50_cases_phenotypes_decisions.py new file mode 100644 index 0000000..76f6c3d --- /dev/null +++ b/api/alembic/versions/9a1c2d3e4f50_cases_phenotypes_decisions.py @@ -0,0 +1,73 @@ +"""cases, phenotypes and triage decisions + +Samples become cases (a proband: a VCF plus observed phenotype terms), and triage gains the two +things it needs: HPO annotations to rank against, and somewhere to record the reviewer's calls. + +Revision ID: 9a1c2d3e4f50 +Revises: 5f2c8e1b9d04 +Create Date: 2026-09-12 09:00:00.000000 + +""" +from collections.abc import Sequence + +import sqlalchemy as sa +from alembic import op + +revision: str = '9a1c2d3e4f50' +down_revision: str | None = '5f2c8e1b9d04' +branch_labels: str | Sequence[str] | None = None +depends_on: str | Sequence[str] | None = None + + +def upgrade() -> None: + op.rename_table('samples', 'cases') + op.alter_column('jobs', 'sample_id', new_column_name='case_id') + + op.create_table( + 'case_phenotypes', + sa.Column('id', sa.Integer(), autoincrement=True, nullable=False), + sa.Column('case_id', sa.UUID(), nullable=False), + sa.Column('hpo_id', sa.String(length=20), nullable=False), + sa.Column('label', sa.String(length=200), nullable=False), + sa.ForeignKeyConstraint(['case_id'], ['cases.id'], ondelete='CASCADE'), + sa.PrimaryKeyConstraint('id'), + sa.UniqueConstraint('case_id', 'hpo_id', name='uq_case_phenotypes_case_term'), + ) + op.create_index(op.f('ix_case_phenotypes_case_id'), 'case_phenotypes', ['case_id']) + + op.create_table( + 'gene_phenotypes', + sa.Column('id', sa.Integer(), autoincrement=True, nullable=False), + sa.Column('gene_symbol', sa.String(length=60), nullable=False), + sa.Column('hpo_id', sa.String(length=20), nullable=False), + sa.Column('hpo_name', sa.String(length=200), nullable=False), + sa.PrimaryKeyConstraint('id'), + sa.UniqueConstraint('gene_symbol', 'hpo_id', name='uq_gene_phenotypes_gene_term'), + ) + op.create_index(op.f('ix_gene_phenotypes_gene_symbol'), 'gene_phenotypes', ['gene_symbol']) + op.create_index(op.f('ix_gene_phenotypes_hpo_id'), 'gene_phenotypes', ['hpo_id']) + + op.create_table( + 'variant_decisions', + sa.Column('id', sa.Integer(), autoincrement=True, nullable=False), + sa.Column('variant_id', sa.Integer(), nullable=False), + sa.Column('state', sa.Enum('shortlisted', 'dismissed', name='decisionstate'), nullable=False), + sa.Column('reason', sa.String(length=120), nullable=True), + sa.Column('note', sa.Text(), nullable=True), + sa.Column('decided_at', sa.DateTime(timezone=True), server_default=sa.text('now()'), nullable=False), + sa.ForeignKeyConstraint(['variant_id'], ['variants.id'], ondelete='CASCADE'), + sa.PrimaryKeyConstraint('id'), + sa.UniqueConstraint('variant_id'), + ) + + +def downgrade() -> None: + op.drop_table('variant_decisions') + sa.Enum(name='decisionstate').drop(op.get_bind(), checkfirst=True) + op.drop_index(op.f('ix_gene_phenotypes_hpo_id'), table_name='gene_phenotypes') + op.drop_index(op.f('ix_gene_phenotypes_gene_symbol'), table_name='gene_phenotypes') + op.drop_table('gene_phenotypes') + op.drop_index(op.f('ix_case_phenotypes_case_id'), table_name='case_phenotypes') + op.drop_table('case_phenotypes') + op.alter_column('jobs', 'case_id', new_column_name='sample_id') + op.rename_table('cases', 'samples') diff --git a/api/app/main.py b/api/app/main.py index 3ec70bb..5edfe07 100644 --- a/api/app/main.py +++ b/api/app/main.py @@ -4,7 +4,7 @@ from fastapi import APIRouter, FastAPI from fastapi.middleware.cors import CORSMiddleware from app.config import settings -from app.routers import jobs, predictions, samples, variants +from app.routers import cases, jobs, phenotypes, variants @asynccontextmanager @@ -23,10 +23,10 @@ app.add_middleware( # The ingress forwards /api/* to this service unchanged, and local dev uses the same prefix. api = APIRouter(prefix="/api") -api.include_router(samples.router, prefix="/samples", tags=["samples"]) +api.include_router(cases.router, prefix="/cases", tags=["cases"]) api.include_router(jobs.router, prefix="/jobs", tags=["jobs"]) api.include_router(variants.router, prefix="/variants", tags=["variants"]) -api.include_router(predictions.router, prefix="/predictions", tags=["predictions"]) +api.include_router(phenotypes.router, prefix="/phenotypes", tags=["phenotypes"]) app.include_router(api) diff --git a/api/app/models.py b/api/app/models.py index 5bed222..ab03607 100644 --- a/api/app/models.py +++ b/api/app/models.py @@ -1,12 +1,23 @@ """SQLAlchemy 2.0 declarative models. -One sample -> many jobs; one job -> many variants; one variant -> one prediction (latest). +A case is a proband: a VCF plus the phenotype terms observed in that patient. One case -> many +jobs (pipeline runs); one job -> many variants; one variant -> one prediction and one decision. """ import enum import uuid from datetime import datetime -from sqlalchemy import DateTime, Enum, Float, ForeignKey, Integer, String, Text, func +from sqlalchemy import ( + DateTime, + Enum, + Float, + ForeignKey, + Integer, + String, + Text, + UniqueConstraint, + func, +) from sqlalchemy.dialects.postgresql import JSONB, UUID from sqlalchemy.orm import DeclarativeBase, Mapped, mapped_column, relationship @@ -22,27 +33,58 @@ class JobStatus(str, enum.Enum): failed = "failed" -class Sample(Base): - __tablename__ = "samples" +class DecisionState(str, enum.Enum): + shortlisted = "shortlisted" + dismissed = "dismissed" + + +class Case(Base): + __tablename__ = "cases" id: Mapped[uuid.UUID] = mapped_column(UUID(as_uuid=True), primary_key=True, default=uuid.uuid4) name: Mapped[str] = mapped_column(String(120), unique=True) vcf_uri: Mapped[str] = mapped_column(Text) assembly: Mapped[str] = mapped_column(String(10), default="GRCh38") created_at: Mapped[datetime] = mapped_column(DateTime(timezone=True), server_default=func.now()) - jobs: Mapped[list["Job"]] = relationship(back_populates="sample") + jobs: Mapped[list["Job"]] = relationship(back_populates="case") + phenotypes: Mapped[list["CasePhenotype"]] = relationship( + back_populates="case", cascade="all, delete-orphan", order_by="CasePhenotype.hpo_id" + ) + + +class CasePhenotype(Base): + """An HPO term observed in this patient; the ranking matches these against gene annotations.""" + + __tablename__ = "case_phenotypes" + __table_args__ = (UniqueConstraint("case_id", "hpo_id", name="uq_case_phenotypes_case_term"),) + id: Mapped[int] = mapped_column(Integer, primary_key=True, autoincrement=True) + case_id: Mapped[uuid.UUID] = mapped_column(ForeignKey("cases.id", ondelete="CASCADE"), index=True) + hpo_id: Mapped[str] = mapped_column(String(20)) + label: Mapped[str] = mapped_column(String(200)) + case: Mapped[Case] = relationship(back_populates="phenotypes") + + +class GenePhenotype(Base): + """HPO's gene-to-phenotype annotations (scripts/load-hpo.py); read-only reference data.""" + + __tablename__ = "gene_phenotypes" + __table_args__ = (UniqueConstraint("gene_symbol", "hpo_id", name="uq_gene_phenotypes_gene_term"),) + id: Mapped[int] = mapped_column(Integer, primary_key=True, autoincrement=True) + gene_symbol: Mapped[str] = mapped_column(String(60), index=True) + hpo_id: Mapped[str] = mapped_column(String(20), index=True) + hpo_name: Mapped[str] = mapped_column(String(200)) class Job(Base): __tablename__ = "jobs" id: Mapped[uuid.UUID] = mapped_column(UUID(as_uuid=True), primary_key=True, default=uuid.uuid4) - sample_id: Mapped[uuid.UUID] = mapped_column(ForeignKey("samples.id", ondelete="CASCADE")) + case_id: Mapped[uuid.UUID] = mapped_column(ForeignKey("cases.id", ondelete="CASCADE")) status: Mapped[JobStatus] = mapped_column(Enum(JobStatus), default=JobStatus.queued) workflow_ref: Mapped[str | None] = mapped_column(String(200)) # Argo workflow name / nf run id vep_version: Mapped[str | None] = mapped_column(String(40)) log: Mapped[str | None] = mapped_column(Text) created_at: Mapped[datetime] = mapped_column(DateTime(timezone=True), server_default=func.now()) finished_at: Mapped[datetime | None] = mapped_column(DateTime(timezone=True)) - sample: Mapped[Sample] = relationship(back_populates="jobs") + case: Mapped[Case] = relationship(back_populates="jobs") variants: Mapped[list["Variant"]] = relationship(back_populates="job") @@ -64,6 +106,7 @@ class Variant(Base): annotations: Mapped[dict] = mapped_column(JSONB, default=dict) # full VEP CSQ record job: Mapped[Job] = relationship(back_populates="variants") prediction: Mapped["Prediction | None"] = relationship(back_populates="variant", uselist=False) + decision: Mapped["VariantDecision | None"] = relationship(back_populates="variant", uselist=False) class Prediction(Base): @@ -75,3 +118,16 @@ class Prediction(Base): score: Mapped[float] = mapped_column(Float) # P(pathogenic) created_at: Mapped[datetime] = mapped_column(DateTime(timezone=True), server_default=func.now()) variant: Mapped[Variant] = relationship(back_populates="prediction") + + +class VariantDecision(Base): + """A reviewer's triage call. No authentication here, so decisions are shared by all visitors.""" + + __tablename__ = "variant_decisions" + id: Mapped[int] = mapped_column(Integer, primary_key=True, autoincrement=True) + variant_id: Mapped[int] = mapped_column(ForeignKey("variants.id", ondelete="CASCADE"), unique=True) + state: Mapped[DecisionState] = mapped_column(Enum(DecisionState)) + reason: Mapped[str | None] = mapped_column(String(120)) + note: Mapped[str | None] = mapped_column(Text) + decided_at: Mapped[datetime] = mapped_column(DateTime(timezone=True), server_default=func.now()) + variant: Mapped[Variant] = relationship(back_populates="decision") diff --git a/api/app/routers/cases.py b/api/app/routers/cases.py new file mode 100644 index 0000000..297e901 --- /dev/null +++ b/api/app/routers/cases.py @@ -0,0 +1,225 @@ +import logging +import uuid +from dataclasses import asdict +from datetime import UTC, datetime + +from fastapi import APIRouter, HTTPException, Query, status +from sqlalchemy import func, select +from sqlalchemy.exc import IntegrityError +from sqlalchemy.orm import selectinload + +from app.db import SessionDep +from app.models import ( + Case, + CasePhenotype, + DecisionState, + Job, + JobStatus, + Variant, + VariantDecision, +) +from app.schemas import ( + CandidateOut, + CandidatePage, + CaseCreate, + CaseOut, + FunnelOut, + JobOut, + ProvenanceOut, + ReportOut, + ScoreOut, +) +from app.services import candidates as case_view +from app.services import events, triage +from app.services.scoring import score_job + +logger = logging.getLogger(__name__) +router = APIRouter() + + +async def _case_or_404(session: SessionDep, case_id: uuid.UUID) -> Case: + case = await case_view.get_case(session, case_id) + if case is None: + raise HTTPException(404, "case not found") + return case + + +async def _shortlisted_counts(session: SessionDep, case_ids: list[uuid.UUID]) -> dict[uuid.UUID, int]: + if not case_ids: + return {} + rows = await session.execute( + select(Job.case_id, func.count(VariantDecision.id)) + .join(Variant, Variant.job_id == Job.id) + .join(VariantDecision, VariantDecision.variant_id == Variant.id) + .where(Job.case_id.in_(case_ids), VariantDecision.state == DecisionState.shortlisted) + .group_by(Job.case_id) + ) + return dict(rows.all()) # type: ignore[arg-type] + + +async def _latest_jobs(session: SessionDep, case_ids: list[uuid.UUID]) -> dict[uuid.UUID, Job]: + if not case_ids: + return {} + jobs = await session.scalars( + select(Job).where(Job.case_id.in_(case_ids)).order_by(Job.created_at.desc()) + ) + latest: dict[uuid.UUID, Job] = {} + for job in jobs: # ordered newest first, so the first one wins + latest.setdefault(job.case_id, job) + return latest + + +def _as_case_out(case: Case, job: Job | None, shortlisted: int) -> CaseOut: + out = CaseOut.model_validate(case) + out.latest_job = JobOut.model_validate(job) if job is not None else None + out.shortlisted = shortlisted + return out + + +@router.get("", response_model=list[CaseOut]) +async def list_cases(session: SessionDep): + cases = ( + await session.scalars( + select(Case).order_by(Case.created_at.desc()).options(selectinload(Case.phenotypes)) + ) + ).all() + ids = [c.id for c in cases] + jobs, counts = await _latest_jobs(session, ids), await _shortlisted_counts(session, ids) + return [_as_case_out(c, jobs.get(c.id), counts.get(c.id, 0)) for c in cases] + + +@router.post("", response_model=CaseOut, status_code=status.HTTP_201_CREATED) +async def create_case(payload: CaseCreate, session: SessionDep): + case = Case( + name=payload.name, + vcf_uri=payload.vcf_uri, + assembly=payload.assembly, + phenotypes=[CasePhenotype(hpo_id=p.hpo_id, label=p.label) for p in payload.phenotypes], + ) + session.add(case) + try: + await session.commit() + except IntegrityError: # cases.name is unique + await session.rollback() + raise HTTPException(409, f"a case named {payload.name!r} already exists") from None + await session.refresh(case, attribute_names=["phenotypes"]) + return _as_case_out(case, None, 0) + + +@router.get("/{case_id}", response_model=CaseOut) +async def get_case(case_id: uuid.UUID, session: SessionDep): + case = await _case_or_404(session, case_id) + counts = await _shortlisted_counts(session, [case_id]) + return _as_case_out(case, await case_view.latest_job(session, case_id), counts.get(case_id, 0)) + + +@router.get("/{case_id}/jobs", response_model=list[JobOut]) +async def list_jobs(case_id: uuid.UUID, session: SessionDep): + result = await session.scalars( + select(Job).where(Job.case_id == case_id).order_by(Job.created_at.desc()) + ) + return result.all() + + +@router.post("/{case_id}/annotate", response_model=JobOut, status_code=status.HTTP_202_ACCEPTED) +async def annotate(case_id: uuid.UUID, session: SessionDep): + case = await _case_or_404(session, case_id) + # Commit `running` before launching: a local run that dies instantly is marked failed by its + # watcher, and a later status write here would overwrite that. + job = Job(case_id=case.id, status=JobStatus.running) + session.add(job) + await session.commit() + try: + job.workflow_ref = await events.launch(job.id, case.vcf_uri, case.assembly) + except events.LaunchError as e: + job.status = JobStatus.failed + job.log = str(e) + job.finished_at = datetime.now(UTC) + await session.commit() + await session.refresh(job) + return job + + +@router.post("/{case_id}/score", response_model=ScoreOut) +async def score(case_id: uuid.UUID, session: SessionDep): + await _case_or_404(session, case_id) + job = await case_view.latest_job(session, case_id, status=JobStatus.succeeded) + if job is None: + raise HTTPException(409, "no finished annotation to score") + try: + scored, version = await score_job(job.id, session) + except Exception as e: + # No registry, no model behind the alias, a model that will not load: all of these are + # the environment being unready, not a bug in the request. Say so rather than throwing 500. + logger.exception("scoring case %s failed", case_id) + raise HTTPException(503, f"could not score with the model: {e}") from e + return ScoreOut(case_id=case_id, scored=scored, model_version=version) + + +@router.get("/{case_id}/candidates", response_model=CandidatePage) +async def list_candidates( + case_id: uuid.UUID, + session: SessionDep, + gene: str | None = None, + impact: str | None = Query(None, pattern="^(HIGH|MODERATE|LOW|MODIFIER)$"), + max_af: float | None = Query(None, ge=0, le=1), + state: str | None = Query(None, pattern="^(shortlisted|dismissed|undecided)$"), + limit: int = Query(50, ge=1, le=500), + offset: int = Query(0, ge=0), +) -> CandidatePage: + case = await _case_or_404(session, case_id) + view = await case_view.build(session, case) + + items = view.candidates + if gene: + items = [c for c in items if (c.variant.gene or "").upper() == gene.upper()] + if impact: + items = [c for c in items if c.variant.impact == impact] + if max_af is not None: + items = [c for c in items if c.variant.gnomad_af is None or c.variant.gnomad_af <= max_af] + if state == "undecided": + items = [c for c in items if c.variant.decision is None] + elif state is not None: + items = [c for c in items if c.variant.decision and c.variant.decision.state.value == state] + + labels = view.labels + return CandidatePage( + # The funnel describes the whole case, not the filtered view. + funnel=FunnelOut(**asdict(view.funnel)), + weights=triage.WEIGHTS, + items=[CandidateOut.from_candidate(c, labels) for c in items[offset : offset + limit]], + total=len(items), + limit=limit, + offset=offset, + ) + + +@router.get("/{case_id}/report", response_model=ReportOut) +async def report(case_id: uuid.UUID, session: SessionDep) -> ReportOut: + case = await _case_or_404(session, case_id) + view = await case_view.build(session, case) + labels = view.labels + + def decided(state: DecisionState) -> list[CandidateOut]: + return [ + CandidateOut.from_candidate(c, labels) + for c in view.candidates + if c.variant.decision is not None and c.variant.decision.state == state + ] + + shortlisted, dismissed = decided(DecisionState.shortlisted), decided(DecisionState.dismissed) + prediction = next((c.variant.prediction for c in view.candidates if c.variant.prediction), None) + return ReportOut( + case=_as_case_out(case, view.job, len(shortlisted)), + funnel=FunnelOut(**asdict(view.funnel)), + generated_at=datetime.now(UTC), + provenance=ProvenanceOut( + job_id=view.job.id if view.job else None, + vep_version=view.job.vep_version if view.job else None, + finished_at=view.job.finished_at if view.job else None, + model_name=prediction.model_name if prediction else None, + model_version=prediction.model_version if prediction else None, + ), + shortlisted=shortlisted, + dismissed=dismissed, + ) diff --git a/api/app/routers/phenotypes.py b/api/app/routers/phenotypes.py new file mode 100644 index 0000000..d37601b --- /dev/null +++ b/api/app/routers/phenotypes.py @@ -0,0 +1,25 @@ +from fastapi import APIRouter, Query +from sqlalchemy import select + +from app.db import SessionDep +from app.models import GenePhenotype +from app.schemas import PhenotypeTerm + +router = APIRouter() + + +@router.get("", response_model=list[PhenotypeTerm]) +async def search( + session: SessionDep, + q: str = Query(min_length=2, max_length=100), + limit: int = Query(20, ge=1, le=100), +) -> list[PhenotypeTerm]: + """Type-ahead over the HPO terms that are annotated to at least one gene.""" + rows = await session.execute( + select(GenePhenotype.hpo_id, GenePhenotype.hpo_name) + .where(GenePhenotype.hpo_name.ilike(f"%{q}%")) + .distinct() + .order_by(GenePhenotype.hpo_name) + .limit(limit) + ) + return [PhenotypeTerm(hpo_id=hpo_id, label=name) for hpo_id, name in rows] diff --git a/api/app/routers/predictions.py b/api/app/routers/predictions.py deleted file mode 100644 index 3732e10..0000000 --- a/api/app/routers/predictions.py +++ /dev/null @@ -1,22 +0,0 @@ -import uuid - -from fastapi import APIRouter, HTTPException - -from app.db import SessionDep -from app.models import Job, JobStatus -from app.schemas import ScoreOut -from app.services.scoring import score_job - -router = APIRouter() - - -@router.post("/score/{job_id}", response_model=ScoreOut) -async def score(job_id: uuid.UUID, session: SessionDep) -> ScoreOut: - """Score every variant of a finished job with the model behind the registry alias.""" - job = await session.get(Job, job_id) - if job is None: - raise HTTPException(404, "job not found") - if job.status != JobStatus.succeeded: - raise HTTPException(409, f"job is {job.status.value}; only succeeded jobs can be scored") - n, version = await score_job(job_id, session) - return ScoreOut(job_id=job_id, scored=n, model_version=version) diff --git a/api/app/routers/samples.py b/api/app/routers/samples.py deleted file mode 100644 index c4b774e..0000000 --- a/api/app/routers/samples.py +++ /dev/null @@ -1,61 +0,0 @@ -import uuid -from datetime import UTC, datetime - -from fastapi import APIRouter, HTTPException, status -from sqlalchemy import select -from sqlalchemy.exc import IntegrityError - -from app.db import SessionDep -from app.models import Job, JobStatus, Sample -from app.schemas import JobOut, SampleCreate, SampleOut -from app.services import events - -router = APIRouter() - - -@router.get("", response_model=list[SampleOut]) -async def list_samples(session: SessionDep): - result = await session.scalars(select(Sample).order_by(Sample.created_at.desc())) - return result.all() - - -@router.post("", response_model=SampleOut, status_code=status.HTTP_201_CREATED) -async def create_sample(payload: SampleCreate, session: SessionDep): - sample = Sample(**payload.model_dump()) - session.add(sample) - try: - await session.commit() - except IntegrityError: # samples.name is unique - await session.rollback() - raise HTTPException(409, f"a sample named {payload.name!r} already exists") from None - await session.refresh(sample) - return sample - - -@router.get("/{sample_id}/jobs", response_model=list[JobOut]) -async def list_jobs(sample_id: uuid.UUID, session: SessionDep): - result = await session.scalars( - select(Job).where(Job.sample_id == sample_id).order_by(Job.created_at.desc()) - ) - return result.all() - - -@router.post("/{sample_id}/annotate", response_model=JobOut, status_code=status.HTTP_202_ACCEPTED) -async def annotate(sample_id: uuid.UUID, session: SessionDep): - sample = await session.get(Sample, sample_id) - if sample is None: - raise HTTPException(404, "sample not found") - # Commit `running` before launching: a local run that dies instantly is marked failed by its - # watcher, and a later status write here would overwrite that. - job = Job(sample_id=sample.id, status=JobStatus.running) - session.add(job) - await session.commit() - try: - job.workflow_ref = await events.launch(job.id, sample.vcf_uri, sample.assembly) - except events.LaunchError as e: - job.status = JobStatus.failed - job.log = str(e) - job.finished_at = datetime.now(UTC) - await session.commit() - await session.refresh(job) - return job diff --git a/api/app/routers/variants.py b/api/app/routers/variants.py index 94f6854..9e88430 100644 --- a/api/app/routers/variants.py +++ b/api/app/routers/variants.py @@ -1,60 +1,63 @@ -import uuid - -from fastapi import APIRouter, Query -from sqlalchemy import Integer, case, cast, func, select +from fastapi import APIRouter, HTTPException +from sqlalchemy import func, select +from sqlalchemy.dialects.postgresql import insert from sqlalchemy.orm import selectinload from app.db import SessionDep -from app.models import Prediction, Variant -from app.schemas import VariantOut, VariantPage +from app.models import Variant, VariantDecision +from app.schemas import CandidateOut, DecisionIn, DecisionOut, VariantDetailOut +from app.services import candidates as case_view +from app.services import triage router = APIRouter() -# Karyotype order (1..22, X, Y, MT) instead of text order, where "10" sorts before "2". -_chrom = func.regexp_replace(Variant.chrom, "^chr", "", "i") -CHROM_ORDER = case( - (_chrom.regexp_match("^[0-9]+$"), cast(_chrom, Integer)), - (_chrom == "X", 23), - (_chrom == "Y", 24), - (_chrom.in_(["M", "MT"]), 25), - else_=26, -) - -@router.get("", response_model=VariantPage) -async def list_variants( - job_id: uuid.UUID, - session: SessionDep, - gene: str | None = None, - impact: str | None = Query(None, pattern="^(HIGH|MODERATE|LOW|MODIFIER)$"), - max_af: float | None = Query(None, ge=0, le=1), - min_score: float | None = Query(None, ge=0, le=1), - limit: int = Query(50, ge=1, le=500), - offset: int = Query(0, ge=0), -) -> VariantPage: - stmt = select(Variant).where(Variant.job_id == job_id) - if gene: - stmt = stmt.where(Variant.gene == gene.upper()) - if impact: - stmt = stmt.where(Variant.impact == impact) - if max_af is not None: - stmt = stmt.where((Variant.gnomad_af.is_(None)) | (Variant.gnomad_af <= max_af)) - if min_score is not None: - stmt = stmt.join(Prediction, Prediction.variant_id == Variant.id).where( - Prediction.score >= min_score +@router.get("/{variant_id}", response_model=VariantDetailOut) +async def get_variant(variant_id: int, session: SessionDep) -> VariantDetailOut: + """Every piece of evidence for one variant, including the raw VEP record.""" + variant = await session.scalar( + select(Variant) + .where(Variant.id == variant_id) + .options( + selectinload(Variant.prediction), + selectinload(Variant.decision), + selectinload(Variant.job), ) + ) + if variant is None: + raise HTTPException(404, "variant not found") + case = await case_view.get_case(session, variant.job.case_id) + terms = [p.hpo_id for p in case.phenotypes] if case else [] + gene_terms = await case_view.gene_terms_for(session, {variant.gene} if variant.gene else set()) - total = await session.scalar(select(func.count()).select_from(stmt.subquery())) - rows = await session.scalars( - stmt.options(selectinload(Variant.prediction)) - # id breaks ties between split multiallelics at one position, keeping pages stable. - .order_by(CHROM_ORDER, Variant.chrom, Variant.pos, Variant.id) - .limit(limit) - .offset(offset) + # evaluate, not rank: the panel must work for a variant that did not make the candidate list. + scored = triage.evaluate(variant, terms, gene_terms) + labels = {p.hpo_id: p.label for p in case.phenotypes} if case else {} + base = CandidateOut.from_candidate(scored, labels).model_dump() + return VariantDetailOut(**base, annotations=variant.annotations or {}) + + +@router.post("/{variant_id}/decision", response_model=DecisionOut) +async def decide(variant_id: int, payload: DecisionIn, session: SessionDep) -> DecisionOut: + if await session.get(Variant, variant_id) is None: + raise HTTPException(404, "variant not found") + stmt = insert(VariantDecision).values( + variant_id=variant_id, state=payload.state, reason=payload.reason, note=payload.note ) - return VariantPage( - items=[VariantOut.model_validate(v) for v in rows], - total=total or 0, - limit=limit, - offset=offset, + # Changing your mind replaces the call rather than failing on the unique constraint. + await session.execute( + stmt.on_conflict_do_update( + index_elements=[VariantDecision.variant_id], + set_={ + "state": stmt.excluded.state, + "reason": stmt.excluded.reason, + "note": stmt.excluded.note, + "decided_at": func.now(), + }, + ) ) + await session.commit() + decision = await session.scalar( + select(VariantDecision).where(VariantDecision.variant_id == variant_id) + ) + return DecisionOut.model_validate(decision) diff --git a/api/app/schemas.py b/api/app/schemas.py index dfc3896..2ddd9f8 100644 --- a/api/app/schemas.py +++ b/api/app/schemas.py @@ -2,12 +2,15 @@ import re import uuid from datetime import datetime from pathlib import PurePosixPath -from typing import Literal +from typing import TYPE_CHECKING, Literal from pydantic import BaseModel, ConfigDict, Field, field_validator from app.config import settings -from app.models import JobStatus +from app.models import DecisionState, JobStatus + +if TYPE_CHECKING: + from app.services.triage import Candidate class ORMModel(BaseModel): @@ -19,11 +22,18 @@ VCF_SUFFIXES = (".vcf", ".vcf.gz", ".vcf.bgz", ".bcf") GCS_URI = re.compile(r"gs://[a-z0-9][a-z0-9._-]{1,220}[a-z0-9]/\S+") -class SampleCreate(BaseModel): +class PhenotypeTerm(ORMModel): + """An HPO term: the id is what ranking matches on, the label is for people.""" + + hpo_id: str = Field(pattern=r"^HP:\d{7}$") + label: str = Field(min_length=1, max_length=200) + + +class CaseCreate(BaseModel): name: str = Field(min_length=1, max_length=120) vcf_uri: str - # Passed to VEP --assembly; the VEP cache must contain it. assembly: Assembly = "GRCh38" + phenotypes: list[PhenotypeTerm] = Field(default_factory=list, max_length=100) @field_validator("vcf_uri") @classmethod @@ -44,17 +54,9 @@ class SampleCreate(BaseModel): return v -class SampleOut(ORMModel): - id: uuid.UUID - name: str - vcf_uri: str - assembly: str - created_at: datetime - - class JobOut(ORMModel): id: uuid.UUID - sample_id: uuid.UUID + case_id: uuid.UUID status: JobStatus workflow_ref: str | None vep_version: str | None @@ -63,6 +65,17 @@ class JobOut(ORMModel): finished_at: datetime | None +class CaseOut(ORMModel): + id: uuid.UUID + name: str + vcf_uri: str + assembly: str + created_at: datetime + phenotypes: list[PhenotypeTerm] = Field(default_factory=list) + latest_job: JobOut | None = None + shortlisted: int = 0 + + class PredictionOut(ORMModel): model_name: str model_version: str @@ -85,14 +98,81 @@ class VariantOut(ORMModel): prediction: PredictionOut | None = None -class ScoreOut(BaseModel): - job_id: uuid.UUID - scored: int - model_version: str +class DecisionIn(BaseModel): + state: DecisionState + reason: str | None = Field(None, max_length=120) + note: str | None = Field(None, max_length=2000) -class VariantPage(BaseModel): - items: list[VariantOut] +class DecisionOut(ORMModel): + state: DecisionState + reason: str | None + note: str | None + decided_at: datetime + + +class CandidateOut(BaseModel): + variant: VariantOut + score: float + components: dict[str, float] + matched_terms: list[PhenotypeTerm] + scored: bool + decision: DecisionOut | None = None + + @classmethod + def from_candidate(cls, candidate: "Candidate", labels: dict[str, str]) -> "CandidateOut": + variant = candidate.variant + return cls( + variant=VariantOut.model_validate(variant), + score=round(candidate.score, 4), + components={name: round(v, 4) for name, v in candidate.components.items()}, + matched_terms=[ + PhenotypeTerm(hpo_id=term, label=labels.get(term, term)) + for term in candidate.matched_terms + ], + scored=candidate.scored, + decision=DecisionOut.model_validate(variant.decision) if variant.decision else None, + ) + + +class VariantDetailOut(CandidateOut): + annotations: dict + + +class FunnelOut(BaseModel): + total: int + rare: int + candidates: int + phenotype_matched: int + + +class CandidatePage(BaseModel): + funnel: FunnelOut + weights: dict[str, float] + items: list[CandidateOut] total: int limit: int offset: int + + +class ProvenanceOut(BaseModel): + job_id: uuid.UUID | None = None + vep_version: str | None = None + finished_at: datetime | None = None + model_name: str | None = None + model_version: str | None = None + + +class ReportOut(BaseModel): + case: CaseOut + funnel: FunnelOut + generated_at: datetime + provenance: ProvenanceOut + shortlisted: list[CandidateOut] + dismissed: list[CandidateOut] + + +class ScoreOut(BaseModel): + case_id: uuid.UUID + scored: int + model_version: str diff --git a/api/app/services/candidates.py b/api/app/services/candidates.py new file mode 100644 index 0000000..4c0272d --- /dev/null +++ b/api/app/services/candidates.py @@ -0,0 +1,79 @@ +"""Assemble a case's ranked candidates: load the latest results, rank them, attach decisions. + +Everything for one job is loaded at once, which is fine for a gene panel or a chromosome — the +size of case this demo handles. A whole genome would need the narrowing pushed into SQL. +""" +import uuid +from dataclasses import dataclass + +from sqlalchemy import select +from sqlalchemy.ext.asyncio import AsyncSession +from sqlalchemy.orm import selectinload + +from app.models import Case, GenePhenotype, Job, JobStatus, Variant +from app.services import triage + +EMPTY_FUNNEL = triage.Funnel(total=0, rare=0, candidates=0, phenotype_matched=0) + + +@dataclass(frozen=True) +class CaseView: + case: Case + job: Job | None + funnel: triage.Funnel + candidates: list[triage.Candidate] + + @property + def labels(self) -> dict[str, str]: + return {p.hpo_id: p.label for p in self.case.phenotypes} + + +async def get_case(session: AsyncSession, case_id: uuid.UUID) -> Case | None: + return await session.scalar( + select(Case).where(Case.id == case_id).options(selectinload(Case.phenotypes)) + ) + + +async def latest_job( + session: AsyncSession, case_id: uuid.UUID, *, status: JobStatus | None = None +) -> Job | None: + stmt = select(Job).where(Job.case_id == case_id).order_by(Job.created_at.desc()).limit(1) + if status is not None: + stmt = stmt.where(Job.status == status) + return await session.scalar(stmt) + + +async def gene_terms_for(session: AsyncSession, genes: set[str]) -> dict[str, set[str]]: + """gene symbol -> the HPO terms annotated to it.""" + if not genes: + return {} + rows = await session.execute( + select(GenePhenotype.gene_symbol, GenePhenotype.hpo_id).where( + GenePhenotype.gene_symbol.in_(genes) + ) + ) + index: dict[str, set[str]] = {} + for gene, hpo_id in rows: + index.setdefault(gene, set()).add(hpo_id) + return index + + +async def build(session: AsyncSession, case: Case) -> CaseView: + job = await latest_job(session, case.id, status=JobStatus.succeeded) + if job is None: + return CaseView(case=case, job=None, funnel=EMPTY_FUNNEL, candidates=[]) + variants = ( + await session.scalars( + select(Variant) + .where(Variant.job_id == job.id) + .options(selectinload(Variant.prediction), selectinload(Variant.decision)) + ) + ).all() + case_terms = [p.hpo_id for p in case.phenotypes] + gene_terms = await gene_terms_for(session, {v.gene for v in variants if v.gene}) + return CaseView( + case=case, + job=job, + funnel=triage.funnel(variants, case_terms, gene_terms), + candidates=triage.rank(variants, case_terms, gene_terms), + ) diff --git a/api/app/services/triage.py b/api/app/services/triage.py new file mode 100644 index 0000000..4ebda62 --- /dev/null +++ b/api/app/services/triage.py @@ -0,0 +1,106 @@ +"""Narrow a case's variants the way a clinical scientist does, and say why. + +The rank is a weighted sum of four parts a reviewer can audit. ClinVar is deliberately not one of +them: it is shown beside the result as independent confirmation, so a variant never ranks highly +merely because ClinVar already called it pathogenic. + +Rarity and consequence *filter* (the usual first pass); phenotype only *ranks*, because a real +diagnosis can sit in a gene nobody has annotated yet and filtering on it would hide exactly that. +""" +from collections.abc import Mapping, Sequence +from dataclasses import dataclass + +from app.models import Variant + +WEIGHTS = {"phenotype": 0.35, "rarity": 0.25, "consequence": 0.20, "model": 0.20} + +RARE_AF = 0.001 +CANDIDATE_IMPACTS = frozenset({"HIGH", "MODERATE"}) +IMPACT_SEVERITY = {"HIGH": 1.0, "MODERATE": 0.6, "LOW": 0.2, "MODIFIER": 0.0} +# Allele frequency ceiling -> score, rarest first. +RARITY_STEPS = ((0.0, 1.0), (0.0001, 0.8), (0.001, 0.5), (0.01, 0.2)) + + +@dataclass(frozen=True) +class Funnel: + """How many variants survive each narrowing step; the headline of the case page.""" + + total: int + rare: int + candidates: int + phenotype_matched: int + + +@dataclass(frozen=True) +class Candidate: + variant: Variant + score: float + components: dict[str, float] + matched_terms: list[str] + scored: bool + + +def rarity_score(af: float | None) -> float: + if af is None: # absent from gnomAD + return 1.0 + for ceiling, score in RARITY_STEPS: + if af <= ceiling: + return score + return 0.0 + + +def consequence_score(impact: str | None) -> float: + return IMPACT_SEVERITY.get(impact or "", 0.0) + + +def phenotype_score( + gene: str | None, case_terms: Sequence[str], gene_terms: Mapping[str, set[str]] +) -> tuple[float, list[str]]: + """What fraction of the patient's terms HPO associates with this gene, and which ones.""" + if not gene or not case_terms: + return 0.0, [] + annotated = gene_terms.get(gene, set()) + matched = [term for term in case_terms if term in annotated] + return len(matched) / len(case_terms), matched + + +def is_rare(variant: Variant) -> bool: + return variant.gnomad_af is None or variant.gnomad_af < RARE_AF + + +def is_candidate(variant: Variant) -> bool: + return is_rare(variant) and variant.impact in CANDIDATE_IMPACTS + + +def funnel( + variants: Sequence[Variant], case_terms: Sequence[str], gene_terms: Mapping[str, set[str]] +) -> Funnel: + rare = [v for v in variants if is_rare(v)] + candidates = [v for v in rare if v.impact in CANDIDATE_IMPACTS] + matched = sum(1 for v in candidates if phenotype_score(v.gene, case_terms, gene_terms)[1]) + return Funnel(len(variants), len(rare), len(candidates), matched) + + +def evaluate( + variant: Variant, case_terms: Sequence[str], gene_terms: Mapping[str, set[str]] +) -> Candidate: + """Score one variant, whether or not it survived the filters.""" + phenotype, matched = phenotype_score(variant.gene, case_terms, gene_terms) + prediction = variant.prediction + components = { + "phenotype": phenotype, + "rarity": rarity_score(variant.gnomad_af), + "consequence": consequence_score(variant.impact), + "model": float(prediction.score) if prediction is not None else 0.0, + } + score = sum(WEIGHTS[name] * value for name, value in components.items()) + return Candidate(variant, score, components, matched, prediction is not None) + + +def rank( + variants: Sequence[Variant], case_terms: Sequence[str], gene_terms: Mapping[str, set[str]] +) -> list[Candidate]: + candidates = [evaluate(v, case_terms, gene_terms) for v in variants if is_candidate(v)] + # id breaks ties, so equal scores do not shuffle between requests. + candidates.sort(key=lambda c: (-c.score, c.variant.id)) + return candidates diff --git a/api/tests/conftest.py b/api/tests/conftest.py index 80692a5..a2adc70 100644 --- a/api/tests/conftest.py +++ b/api/tests/conftest.py @@ -70,7 +70,8 @@ async def db(migrated_db: None) -> AsyncIterator[None]: async with engine.begin() as conn: await conn.execute( - text("TRUNCATE samples, jobs, variants, predictions RESTART IDENTITY CASCADE") + text("TRUNCATE cases, case_phenotypes, gene_phenotypes, jobs, variants, " + "predictions, variant_decisions RESTART IDENTITY CASCADE") ) yield diff --git a/api/tests/factories.py b/api/tests/factories.py index 8fc481e..43c4928 100644 --- a/api/tests/factories.py +++ b/api/tests/factories.py @@ -2,20 +2,47 @@ import uuid from typing import Any from app.db import SessionLocal -from app.models import Job, JobStatus, Sample, Variant +from app.models import Case, CasePhenotype, GenePhenotype, Job, JobStatus, Prediction, Variant + +VARIANT_DEFAULTS: dict[str, Any] = { + "chrom": "22", "pos": 1, "ref": "A", "alt": "G", "gene": "NF2", + "impact": "HIGH", "consequence": "frameshift_variant", "gnomad_af": None, "annotations": {}, +} -async def seed_job( - variants: list[dict[str, Any]], status: JobStatus = JobStatus.succeeded -) -> uuid.UUID: - """Insert a sample, a job and its variants; each variant dict overrides the defaults.""" +async def seed_case( + *, + phenotypes: list[tuple[str, str]] | None = None, + variants: list[dict[str, Any]] | None = None, + gene_terms: dict[str, list[tuple[str, str]]] | None = None, + status: JobStatus = JobStatus.succeeded, + name: str | None = None, +) -> tuple[uuid.UUID, uuid.UUID]: + """Insert a case, its phenotypes, a job and its variants. Returns (case_id, job_id). + + `variants` entries override VARIANT_DEFAULTS; a "score" key becomes a Prediction. + """ async with SessionLocal() as s: - sample = Sample(name=f"s-{uuid.uuid4()}", vcf_uri="gs://b/x.vcf.gz", assembly="GRCh38") - job = Job(sample=sample, status=status) - rows = [ - Variant(job=job, **{"chrom": "22", "pos": 1, "ref": "A", "alt": "G", "annotations": {}} | v) - for v in variants - ] - s.add_all([sample, job, *rows]) + case = Case( + name=name or f"case-{uuid.uuid4()}", + vcf_uri="gs://bucket/proband.vcf.gz", + assembly="GRCh38", + phenotypes=[CasePhenotype(hpo_id=hpo, label=label) for hpo, label in (phenotypes or [])], + ) + job = Job(case=case, status=status, vep_version="113.0") + s.add_all([case, job]) + for gene, terms in (gene_terms or {}).items(): + s.add_all( + GenePhenotype(gene_symbol=gene, hpo_id=hpo, hpo_name=label) for hpo, label in terms + ) + for spec in variants or []: + fields = VARIANT_DEFAULTS | spec + score = fields.pop("score", None) + variant = Variant(job=job, **fields) + s.add(variant) + if score is not None: + await s.flush() + s.add(Prediction(variant_id=variant.id, model_name="rarelens-pathogenicity", + model_version="demo", score=float(score))) await s.commit() - return job.id + return case.id, job.id diff --git a/api/tests/test_annotate.py b/api/tests/test_annotate.py index 0ca2ac2..d6b658e 100644 --- a/api/tests/test_annotate.py +++ b/api/tests/test_annotate.py @@ -8,9 +8,9 @@ from app.config import settings from app.services import events -async def new_sample(client: AsyncClient) -> str: +async def new_case(client: AsyncClient) -> str: r = await client.post( - "/api/samples", + "/api/cases", json={"name": f"s-{uuid.uuid4()}", "vcf_uri": "gs://bucket/x.vcf.gz", "assembly": "GRCh37"}, ) assert r.status_code == 201, r.text @@ -23,9 +23,9 @@ async def test_local_without_nextflow_fails_fast_with_instructions( ) -> None: monkeypatch.setattr(settings, "pubsub_topic", None) monkeypatch.setattr(events.shutil, "which", lambda _: None) - sample_id = await new_sample(client) + case_id = await new_case(client) - r = await client.post(f"/api/samples/{sample_id}/annotate") + r = await client.post(f"/api/cases/{case_id}/annotate") assert r.status_code == 202 job = r.json() assert job["status"] == "failed" @@ -71,9 +71,9 @@ async def test_local_run_is_watched_and_a_crash_marks_the_job_failed( return FakeProcess(1, b"ERROR ~ VEP cache not found") monkeypatch.setattr(events.asyncio, "create_subprocess_exec", fake_exec) - sample_id = await new_sample(client) + case_id = await new_case(client) - r = await client.post(f"/api/samples/{sample_id}/annotate") + r = await client.post(f"/api/cases/{case_id}/annotate") assert r.status_code == 202 await events.drain() @@ -117,9 +117,9 @@ async def test_pubsub_publishes_to_the_full_topic_path( monkeypatch.setattr(settings, "pubsub_topic", "vcf-uploaded") monkeypatch.setattr(settings, "gcp_project", "my-proj") monkeypatch.setattr(events, "_publisher", lambda: publisher) - sample_id = await new_sample(client) + case_id = await new_case(client) - r = await client.post(f"/api/samples/{sample_id}/annotate") + r = await client.post(f"/api/cases/{case_id}/annotate") job = r.json() assert (job["status"], job["workflow_ref"]) == ("running", "pubsub:msg-123") [(topic, data)] = publisher.published @@ -134,9 +134,9 @@ async def test_pubsub_failure_marks_the_job_failed( monkeypatch.setattr(settings, "pubsub_topic", "vcf-uploaded") monkeypatch.setattr(settings, "gcp_project", "my-proj") monkeypatch.setattr(events, "_publisher", lambda: FakePublisher(RuntimeError("403 denied"))) - sample_id = await new_sample(client) + case_id = await new_case(client) - job = (await client.post(f"/api/samples/{sample_id}/annotate")).json() + job = (await client.post(f"/api/cases/{case_id}/annotate")).json() assert job["status"] == "failed" assert "403 denied" in job["log"] @@ -159,9 +159,9 @@ async def test_the_pipeline_gets_its_own_database_url( return FakeProcess(0, b"") monkeypatch.setattr(events.asyncio, "create_subprocess_exec", fake_exec) - sample_id = await new_sample(client) + case_id = await new_case(client) - await client.post(f"/api/samples/{sample_id}/annotate") + await client.post(f"/api/cases/{case_id}/annotate") await events.drain() assert launched["env"]["DATABASE_URL"] == "postgresql+asyncpg://u:p@host.docker.internal:5432/db" @@ -183,9 +183,9 @@ async def test_progress_is_recorded_while_the_pipeline_runs( ]) monkeypatch.setattr(events.asyncio, "create_subprocess_exec", fake_exec) - sample_id = await new_sample(client) + case_id = await new_case(client) - r = await client.post(f"/api/samples/{sample_id}/annotate") + r = await client.post(f"/api/cases/{case_id}/annotate") await events.drain() job = (await client.get(f"/api/jobs/{r.json()['id']}")).json() @@ -196,12 +196,12 @@ async def test_progress_is_recorded_while_the_pipeline_runs( @pytest.mark.usefixtures("db") async def test_progress_never_overwrites_a_finished_job(client: AsyncClient) -> None: from app.db import SessionLocal - from app.models import Job, JobStatus, Sample + from app.models import Case, Job, JobStatus async with SessionLocal() as s: - sample = Sample(name=f"s-{uuid.uuid4()}", vcf_uri="gs://b/x.vcf.gz", assembly="GRCh38") - job = Job(sample=sample, status=JobStatus.succeeded) - s.add_all([sample, job]) + case = Case(name=f"c-{uuid.uuid4()}", vcf_uri="gs://bucket/x.vcf.gz", assembly="GRCh38") + job = Job(case=case, status=JobStatus.succeeded) + s.add_all([case, job]) await s.commit() job_id = job.id diff --git a/api/tests/test_cases.py b/api/tests/test_cases.py new file mode 100644 index 0000000..8e1ab72 --- /dev/null +++ b/api/tests/test_cases.py @@ -0,0 +1,137 @@ +"""The triage workflow: a case with a phenotype, ranked candidates, decisions, a report.""" +import uuid + +import pytest +from factories import seed_case +from httpx import AsyncClient + +NF2_TERMS = [("HP:0000365", "Hearing impairment"), ("HP:0009592", "Vestibular schwannoma")] +CASE_TERMS = [*NF2_TERMS, ("HP:0002321", "Vertigo")] + + +async def a_case_with_candidates() -> tuple[uuid.UUID, uuid.UUID]: + return await seed_case( + phenotypes=CASE_TERMS, + gene_terms={"NF2": NF2_TERMS}, + variants=[ + {"id": 1, "gene": "NF2", "pos": 1000, "impact": "HIGH", "score": 0.94, + "clinvar_sig": "pathogenic"}, + {"id": 2, "gene": "CHEK2", "pos": 2000, "impact": "MODERATE", "gnomad_af": 0.0004, + "score": 0.55}, + {"id": 3, "gene": "TTN", "pos": 3000, "impact": "MODIFIER", "gnomad_af": None}, + {"id": 4, "gene": "APOE", "pos": 4000, "impact": "HIGH", "gnomad_af": 0.3}, + ], + ) + + +@pytest.mark.usefixtures("db") +async def test_a_case_keeps_its_phenotype(client: AsyncClient) -> None: + body = { + "name": "PROBAND-01", + "vcf_uri": "gs://bucket/proband.vcf.gz", + "phenotypes": [{"hpo_id": h, "label": lab} for h, lab in NF2_TERMS], + } + r = await client.post("/api/cases", json=body) + assert r.status_code == 201, r.text + assert [p["label"] for p in r.json()["phenotypes"]] == [lab for _, lab in NF2_TERMS] + + listed = (await client.get("/api/cases")).json() + assert listed[0]["name"] == "PROBAND-01" + assert listed[0]["shortlisted"] == 0 + + +@pytest.mark.usefixtures("db") +async def test_the_funnel_shows_the_narrowing(client: AsyncClient) -> None: + case_id, _ = await a_case_with_candidates() + page = (await client.get(f"/api/cases/{case_id}/candidates")).json() + # 4 variants -> 3 rare -> 2 rare and coding -> 1 of those in a phenotype-matched gene + assert page["funnel"] == {"total": 4, "rare": 3, "candidates": 2, "phenotype_matched": 1} + + +@pytest.mark.usefixtures("db") +async def test_the_phenotype_matched_variant_ranks_first_with_its_reasons( + client: AsyncClient, +) -> None: + case_id, _ = await a_case_with_candidates() + page = (await client.get(f"/api/cases/{case_id}/candidates")).json() + + assert [c["variant"]["gene"] for c in page["items"]] == ["NF2", "CHEK2"] + top = page["items"][0] + assert [t["label"] for t in top["matched_terms"]] == ["Hearing impairment", "Vestibular schwannoma"] + assert top["components"]["phenotype"] == pytest.approx(2 / 3, abs=1e-4) + assert top["components"]["rarity"] == 1.0 + assert page["weights"]["phenotype"] == 0.35 + # ClinVar is evidence, not an input to the rank. + assert top["variant"]["clinvar_sig"] == "pathogenic" + + +@pytest.mark.usefixtures("db") +async def test_a_reviewer_decides_and_the_decision_sticks(client: AsyncClient) -> None: + case_id, _ = await a_case_with_candidates() + variant_id = (await client.get(f"/api/cases/{case_id}/candidates")).json()["items"][0]["variant"]["id"] + + r = await client.post( + f"/api/variants/{variant_id}/decision", + json={"state": "shortlisted", "reason": "fits the phenotype", "note": "confirm by Sanger"}, + ) + assert r.status_code == 200, r.text + + # Changing your mind replaces the decision rather than failing. + r = await client.post(f"/api/variants/{variant_id}/decision", json={"state": "dismissed"}) + assert r.status_code == 200 + detail = (await client.get(f"/api/variants/{variant_id}")).json() + assert detail["decision"]["state"] == "dismissed" + assert (await client.get("/api/cases")).json()[0]["shortlisted"] == 0 + + +@pytest.mark.usefixtures("db") +async def test_the_variant_panel_carries_the_evidence(client: AsyncClient) -> None: + case_id, _ = await a_case_with_candidates() + variant_id = (await client.get(f"/api/cases/{case_id}/candidates")).json()["items"][0]["variant"]["id"] + + detail = (await client.get(f"/api/variants/{variant_id}")).json() + assert detail["variant"]["hgvsp"] is None or isinstance(detail["variant"]["hgvsp"], str) + assert detail["score"] > 0 + assert [t["hpo_id"] for t in detail["matched_terms"]] == [h for h, _ in NF2_TERMS] + assert "annotations" in detail + + +@pytest.mark.usefixtures("db") +async def test_the_report_is_the_decision_trail(client: AsyncClient) -> None: + case_id, _ = await a_case_with_candidates() + items = (await client.get(f"/api/cases/{case_id}/candidates")).json()["items"] + await client.post(f"/api/variants/{items[0]['variant']['id']}/decision", + json={"state": "shortlisted", "reason": "fits the phenotype"}) + await client.post(f"/api/variants/{items[1]['variant']['id']}/decision", + json={"state": "dismissed", "reason": "gene unrelated to phenotype"}) + + report = (await client.get(f"/api/cases/{case_id}/report")).json() + assert report["funnel"]["candidates"] == 2 + assert [v["variant"]["gene"] for v in report["shortlisted"]] == ["NF2"] + assert report["shortlisted"][0]["decision"]["reason"] == "fits the phenotype" + assert [v["variant"]["gene"] for v in report["dismissed"]] == ["CHEK2"] + assert report["provenance"]["vep_version"] == "113.0" + + +@pytest.mark.usefixtures("db") +async def test_candidates_can_still_be_filtered(client: AsyncClient) -> None: + case_id, _ = await a_case_with_candidates() + page = (await client.get(f"/api/cases/{case_id}/candidates", params={"gene": "chek2"})).json() + assert [c["variant"]["gene"] for c in page["items"]] == ["CHEK2"] + # The funnel still describes the whole case, not the filtered view. + assert page["funnel"]["total"] == 4 + + +@pytest.mark.usefixtures("db") +async def test_phenotype_search_backs_the_picker(client: AsyncClient) -> None: + await seed_case(gene_terms={"NF2": NF2_TERMS}, variants=[]) + found = (await client.get("/api/phenotypes", params={"q": "vestibular"})).json() + assert found == [{"hpo_id": "HP:0009592", "label": "Vestibular schwannoma"}] + + +@pytest.mark.usefixtures("db") +async def test_candidates_before_the_pipeline_has_run(client: AsyncClient) -> None: + r = await client.post("/api/cases", json={"name": "empty", "vcf_uri": "gs://bucket/x.vcf.gz"}) + page = (await client.get(f"/api/cases/{r.json()['id']}/candidates")).json() + assert page["items"] == [] + assert page["funnel"]["total"] == 0 diff --git a/api/tests/test_cloudrun.py b/api/tests/test_cloudrun.py index 7309844..e5e6b2f 100644 --- a/api/tests/test_cloudrun.py +++ b/api/tests/test_cloudrun.py @@ -26,9 +26,9 @@ class FakeJobsClient: return self.result -async def new_sample(client: AsyncClient) -> str: +async def new_case(client: AsyncClient) -> str: r = await client.post( - "/api/samples", + "/api/cases", json={"name": f"s-{uuid.uuid4()}", "vcf_uri": "gs://bucket/x.vcf.gz", "assembly": "GRCh37"}, ) assert r.status_code == 201, r.text @@ -49,9 +49,9 @@ async def test_annotate_executes_the_job_with_pipeline_arguments( ) -> None: jobs = FakeJobsClient(FakeOperation("projects/p/locations/l/executions/rarelens-nextflow-abc12")) monkeypatch.setattr(events, "_jobs_client", lambda: jobs) - sample_id = await new_sample(client) + case_id = await new_case(client) - r = await client.post(f"/api/samples/{sample_id}/annotate") + r = await client.post(f"/api/cases/{case_id}/annotate") assert r.status_code == 202 job = r.json() assert job["status"] == "running" @@ -74,9 +74,9 @@ async def test_a_failed_execution_call_marks_the_job_failed( monkeypatch.setattr( events, "_jobs_client", lambda: FakeJobsClient(RuntimeError("403 permission denied")) ) - sample_id = await new_sample(client) + case_id = await new_case(client) - job = (await client.post(f"/api/samples/{sample_id}/annotate")).json() + job = (await client.post(f"/api/cases/{case_id}/annotate")).json() assert job["status"] == "failed" assert "403 permission denied" in job["log"] @@ -91,7 +91,7 @@ async def test_cloud_run_job_takes_precedence_over_pubsub( monkeypatch.setattr( events, "_publisher", lambda: pytest.fail("Pub/Sub must not be used in the serverless track") ) - sample_id = await new_sample(client) + case_id = await new_case(client) - assert (await client.post(f"/api/samples/{sample_id}/annotate")).json()["status"] == "running" + assert (await client.post(f"/api/cases/{case_id}/annotate")).json()["status"] == "running" assert len(jobs.requests) == 1 diff --git a/api/tests/test_routes.py b/api/tests/test_routes.py index b8d133b..7ad38c6 100644 --- a/api/tests/test_routes.py +++ b/api/tests/test_routes.py @@ -5,8 +5,8 @@ from httpx import AsyncClient @pytest.mark.usefixtures("db") async def test_resources_live_under_api_prefix(client: AsyncClient) -> None: # The ingress forwards /api/* unchanged, so the app itself must serve that prefix. - assert (await client.get("/api/samples")).status_code == 200 - assert (await client.get("/samples")).status_code == 404 + assert (await client.get("/api/cases")).status_code == 200 + assert (await client.get("/cases")).status_code == 404 async def test_health_stays_at_root_for_probes(client: AsyncClient) -> None: diff --git a/api/tests/test_samples.py b/api/tests/test_samples.py deleted file mode 100644 index 0d4b8de..0000000 --- a/api/tests/test_samples.py +++ /dev/null @@ -1,18 +0,0 @@ -import pytest -from httpx import AsyncClient - - -@pytest.mark.usefixtures("db") -async def test_duplicate_sample_name_is_409_not_500(client: AsyncClient) -> None: - body = {"name": "HG002", "vcf_uri": "gs://bucket/hg002.vcf.gz"} - assert (await client.post("/api/samples", json=body)).status_code == 201 - r = await client.post("/api/samples", json=body) - assert r.status_code == 409 - assert "HG002" in r.json()["detail"] - - -async def test_unknown_assembly_is_rejected(client: AsyncClient) -> None: - r = await client.post( - "/api/samples", json={"name": "a", "vcf_uri": "gs://bucket/a.vcf.gz", "assembly": "hg19"} - ) - assert r.status_code == 422 diff --git a/api/tests/test_scoring.py b/api/tests/test_scoring.py index cb85d9a..1b22921 100644 --- a/api/tests/test_scoring.py +++ b/api/tests/test_scoring.py @@ -4,11 +4,12 @@ import uuid import numpy as np import pandas as pd import pytest +from factories import seed_case from httpx import AsyncClient from sqlalchemy import select from app.db import SessionLocal -from app.models import Job, JobStatus, Prediction, Sample, Variant +from app.models import JobStatus, Prediction, Variant from app.services import scoring @@ -40,20 +41,16 @@ class FakeModel: return np.full(len(frame), self.score) -async def make_job(status: JobStatus, n_variants: int) -> uuid.UUID: - async with SessionLocal() as s: - sample = Sample(name=f"s-{uuid.uuid4()}", vcf_uri="gs://b/x.vcf.gz", assembly="GRCh38") - job = Job(sample=sample, status=status) - s.add_all([sample, job, *(variant(job=job, pos=i + 1) for i in range(n_variants))]) - await s.commit() - return job.id +async def make_case(status: JobStatus, n_variants: int) -> tuple[uuid.UUID, uuid.UUID]: + return await seed_case( + status=status, + variants=[{"pos": i + 1, "annotations": {}} for i in range(n_variants)], + ) async def predictions(job_id: uuid.UUID) -> list[Prediction]: async with SessionLocal() as s: - rows = await s.scalars( - select(Prediction).join(Variant).where(Variant.job_id == job_id) - ) + rows = await s.scalars(select(Prediction).join(Variant).where(Variant.job_id == job_id)) return list(rows) @@ -61,15 +58,15 @@ async def predictions(job_id: uuid.UUID) -> list[Prediction]: async def test_scoring_twice_updates_instead_of_failing( client: AsyncClient, monkeypatch: pytest.MonkeyPatch ) -> None: - job_id = await make_job(JobStatus.succeeded, n_variants=3) + case_id, job_id = await make_case(JobStatus.succeeded, n_variants=3) monkeypatch.setattr(scoring, "load_model", lambda: (FakeModel(0.9), "7")) - r = await client.post(f"/api/predictions/score/{job_id}") + r = await client.post(f"/api/cases/{case_id}/score") assert r.status_code == 200, r.text - assert r.json() == {"job_id": str(job_id), "scored": 3, "model_version": "7"} + assert r.json() == {"case_id": str(case_id), "scored": 3, "model_version": "7"} monkeypatch.setattr(scoring, "load_model", lambda: (FakeModel(0.2), "8")) - r = await client.post(f"/api/predictions/score/{job_id}") + r = await client.post(f"/api/cases/{case_id}/score") assert r.status_code == 200, r.text preds = await predictions(job_id) @@ -78,13 +75,28 @@ async def test_scoring_twice_updates_instead_of_failing( @pytest.mark.usefixtures("db") -async def test_scoring_unknown_job_is_404(client: AsyncClient) -> None: - r = await client.post(f"/api/predictions/score/{uuid.uuid4()}") +async def test_scoring_an_unknown_case_is_404(client: AsyncClient) -> None: + r = await client.post(f"/api/cases/{uuid.uuid4()}/score") assert r.status_code == 404 @pytest.mark.usefixtures("db") -async def test_scoring_unfinished_job_is_409(client: AsyncClient) -> None: - job_id = await make_job(JobStatus.running, n_variants=1) - r = await client.post(f"/api/predictions/score/{job_id}") +async def test_scoring_before_the_annotation_finishes_is_409(client: AsyncClient) -> None: + case_id, _ = await make_case(JobStatus.running, n_variants=1) + r = await client.post(f"/api/cases/{case_id}/score") assert r.status_code == 409 + + +@pytest.mark.usefixtures("db") +async def test_an_unreachable_model_registry_is_explained_not_a_500( + client: AsyncClient, monkeypatch: pytest.MonkeyPatch +) -> None: + case_id, _ = await make_case(JobStatus.succeeded, n_variants=1) + + def unreachable() -> tuple: + raise ConnectionError("connection refused to http://localhost:5000") + + monkeypatch.setattr(scoring, "load_model", unreachable) + r = await client.post(f"/api/cases/{case_id}/score") + assert r.status_code == 503 + assert "connection refused" in r.json()["detail"] diff --git a/api/tests/test_security.py b/api/tests/test_security.py index 86b6bdd..e122d62 100644 --- a/api/tests/test_security.py +++ b/api/tests/test_security.py @@ -2,7 +2,7 @@ import pytest from httpx import AsyncClient from pydantic import ValidationError -from app.schemas import SampleCreate +from app.schemas import CaseCreate @pytest.mark.parametrize( @@ -16,7 +16,7 @@ from app.schemas import SampleCreate ], ) def test_vcf_uri_accepts_gcs_objects_and_files_under_the_data_root(uri: str) -> None: - assert SampleCreate(name="s", vcf_uri=uri).vcf_uri == uri + assert CaseCreate(name="s", vcf_uri=uri).vcf_uri == uri @pytest.mark.parametrize( @@ -36,11 +36,11 @@ def test_vcf_uri_accepts_gcs_objects_and_files_under_the_data_root(uri: str) -> ) def test_vcf_uri_rejects_everything_else(uri: str) -> None: with pytest.raises(ValidationError): - SampleCreate(name="s", vcf_uri=uri) + CaseCreate(name="s", vcf_uri=uri) async def test_bad_vcf_uri_is_422_at_the_api(client: AsyncClient) -> None: - r = await client.post("/api/samples", json={"name": "s", "vcf_uri": "/etc/passwd"}) + r = await client.post("/api/cases", json={"name": "s", "vcf_uri": "/etc/passwd"}) assert r.status_code == 422 diff --git a/api/tests/test_triage.py b/api/tests/test_triage.py new file mode 100644 index 0000000..8652e58 --- /dev/null +++ b/api/tests/test_triage.py @@ -0,0 +1,101 @@ +"""Ranking is the scientific claim this app makes, so it is tested as pure logic.""" +import pytest + +from app.models import Prediction, Variant +from app.services import triage + + +def variant(**kw: object) -> Variant: + fields: dict = { + "id": 1, "chrom": "22", "pos": 100, "ref": "A", "alt": "G", + "gene": "NF2", "impact": "HIGH", "consequence": "frameshift_variant", + "gnomad_af": None, "clinvar_sig": None, "annotations": {}, + } + fields.update(kw) + score = fields.pop("score", None) + v = Variant(**fields) + if score is not None: + v.prediction = Prediction(model_name="m", model_version="1", score=float(score)) + return v + + +def test_weights_sum_to_one() -> None: + assert sum(triage.WEIGHTS.values()) == pytest.approx(1.0) + + +@pytest.mark.parametrize( + ("af", "expected"), + [(None, 1.0), (0.0, 1.0), (0.00005, 0.8), (0.0005, 0.5), (0.005, 0.2), (0.05, 0.0)], +) +def test_rarity_rewards_absence_from_gnomad(af: float | None, expected: float) -> None: + assert triage.rarity_score(af) == expected + + +@pytest.mark.parametrize( + ("impact", "expected"), + [("HIGH", 1.0), ("MODERATE", 0.6), ("LOW", 0.2), ("MODIFIER", 0.0), (None, 0.0), ("?", 0.0)], +) +def test_consequence_severity(impact: str | None, expected: float) -> None: + assert triage.consequence_score(impact) == expected + + +def test_phenotype_match_is_the_fraction_of_the_patients_terms() -> None: + gene_terms = {"NF2": {"HP:0000365", "HP:0009592"}} + case_terms = ["HP:0000365", "HP:0009592", "HP:0002321", "HP:0000598"] + score, matched = triage.phenotype_score("NF2", case_terms, gene_terms) + assert score == 0.5 + assert matched == ["HP:0000365", "HP:0009592"] + + +def test_phenotype_match_is_zero_for_genes_hpo_has_never_annotated() -> None: + assert triage.phenotype_score("NOVEL1", ["HP:0000365"], {}) == (0.0, []) + + +def test_phenotype_match_is_zero_when_no_phenotype_was_entered() -> None: + assert triage.phenotype_score("NF2", [], {"NF2": {"HP:0000365"}}) == (0.0, []) + + +def test_the_funnel_counts_each_narrowing_step() -> None: + variants = [ + variant(id=1, gnomad_af=None, impact="HIGH", gene="NF2"), # rare, coding, matched + variant(id=2, gnomad_af=0.0002, impact="MODERATE", gene="CHEK2"), # rare, coding + variant(id=3, gnomad_af=0.3, impact="HIGH", gene="NF2"), # common + variant(id=4, gnomad_af=None, impact="MODIFIER", gene="NF2"), # rare, non-coding + ] + funnel = triage.funnel(variants, case_terms=["HP:0000365"], gene_terms={"NF2": {"HP:0000365"}}) + assert (funnel.total, funnel.rare, funnel.candidates, funnel.phenotype_matched) == (4, 3, 2, 1) + + +def test_the_diagnosis_outranks_the_noise() -> None: + gene_terms = {"NF2": {"HP:0000365", "HP:0009592"}} + case_terms = ["HP:0000365", "HP:0009592"] + diagnosis = variant(id=1, gene="NF2", impact="HIGH", gnomad_af=None, score=0.94) + plausible = variant(id=2, gene="CHEK2", impact="MODERATE", gnomad_af=0.0004, score=0.55) + noise = variant(id=3, gene="TTN", impact="MODERATE", gnomad_af=0.0009, score=0.10) + + ranked = triage.rank([noise, plausible, diagnosis], case_terms, gene_terms) + assert [c.variant.id for c in ranked] == [1, 2, 3] + top = ranked[0] + assert top.matched_terms == case_terms + assert top.components["phenotype"] == 1.0 + assert top.score == pytest.approx(0.35 + 0.25 + 0.20 + 0.20 * 0.94) + + +def test_an_unscored_variant_still_ranks_and_says_so() -> None: + [candidate] = triage.rank([variant(id=1, gnomad_af=None)], [], {}) + assert candidate.components["model"] == 0.0 + assert candidate.scored is False + + +def test_common_and_non_coding_variants_are_not_candidates() -> None: + variants = [ + variant(id=1, gnomad_af=0.2, impact="HIGH"), + variant(id=2, gnomad_af=None, impact="MODIFIER"), + ] + assert triage.rank(variants, [], {}) == [] + + +def test_ranking_is_deterministic_for_equal_scores() -> None: + a = variant(id=7, gene="AAA", chrom="1", pos=10, gnomad_af=None) + b = variant(id=3, gene="BBB", chrom="1", pos=10, gnomad_af=None) + assert [c.variant.id for c in triage.rank([a, b], [], {})] == [3, 7] diff --git a/api/tests/test_variants.py b/api/tests/test_variants.py deleted file mode 100644 index 386cadd..0000000 --- a/api/tests/test_variants.py +++ /dev/null @@ -1,46 +0,0 @@ -import uuid - -import pytest -from factories import seed_job -from httpx import AsyncClient - - -@pytest.mark.parametrize( - "paging", [{"limit": 0}, {"limit": -1}, {"limit": 501}, {"offset": -1}] -) -async def test_out_of_range_paging_is_422_not_500(client: AsyncClient, paging: dict) -> None: - r = await client.get("/api/variants", params={"job_id": str(uuid.uuid4()), **paging}) - assert r.status_code == 422 - - -async def positions(client: AsyncClient, job_id: uuid.UUID) -> list[tuple[str, int, str]]: - r = await client.get("/api/variants", params={"job_id": str(job_id)}) - assert r.status_code == 200, r.text - return [(v["chrom"], v["pos"], v["alt"]) for v in r.json()["items"]] - - -@pytest.mark.usefixtures("db") -async def test_chromosomes_sort_naturally(client: AsyncClient) -> None: - job_id = await seed_job([{"chrom": c} for c in ["10", "MT", "2", "X", "chr3", "1", "Y"]]) - assert [c for c, _, _ in await positions(client, job_id)] == [ - "1", "2", "chr3", "10", "X", "Y", "MT", - ] - - -@pytest.mark.usefixtures("db") -async def test_same_position_keeps_insertion_order_across_pages(client: AsyncClient) -> None: - # Split multiallelics share chrom/pos; without a tiebreak, pages can repeat or skip rows. - job_id = await seed_job([{"pos": 5, "alt": a} for a in "CGT"]) - assert [a for _, _, a in await positions(client, job_id)] == ["C", "G", "T"] - - -@pytest.mark.usefixtures("db") -async def test_long_vep_strings_are_stored(client: AsyncClient) -> None: - clin_sig = ",".join(["conflicting_classifications_of_pathogenicity"] * 8) - consequence = ( - "splice_region_variant&splice_polypyrimidine_tract_variant&intron_variant" - "&NMD_transcript_variant&non_coding_transcript_variant" - ) - job_id = await seed_job([{"clinvar_sig": clin_sig, "consequence": consequence}]) - [v] = (await client.get("/api/variants", params={"job_id": str(job_id)})).json()["items"] - assert (v["clinvar_sig"], v["consequence"]) == (clin_sig, consequence) diff --git a/data/README.md b/data/README.md index c062cec..e3ed81a 100644 --- a/data/README.md +++ b/data/README.md @@ -1,6 +1,6 @@ # Test data -`make data` fetches the demo slice automatically. Provenance, licences, citations and the +`make data` fetches the demo slice and `make demo-case` builds the simulated proband. Provenance, licences, citations and the evaluation plan live in [../docs/data.md](../docs/data.md). No patient data. Use public sources only: diff --git a/docs/architecture.md b/docs/architecture.md index c807654..0299027 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -2,8 +2,9 @@ ```mermaid flowchart LR - U[Scientist] -->|browser| W[SvelteKit web] + U[Scientist] -->|phenotype + VCF| W[SvelteKit web] W -->|REST /api| A[FastAPI] + H[(HPO gene-phenotype annotations)] --> A A --> P[(PostgreSQL / Cloud SQL)] A -->|publish vcf-uploaded| Q[Pub/Sub] Q --> E[Argo Events sensor] @@ -12,6 +13,7 @@ flowchart LR B -->|reads VCF, VEP cache| G[(GCS bucket)] B -->|writes variants, marks job succeeded| P AW -.->|exit handler marks job failed| P + A -->|rank: phenotype, rarity, consequence, model| C[Ranked candidates -> decisions -> report] A -->|models:/rarelens-pathogenicity@production| M[MLflow registry] T[ml/train.py] --> M GH[GitHub Actions] -->|images via WIF| AR[Artifact Registry] @@ -19,6 +21,33 @@ flowchart LR R --> CD[ArgoCD] --> K[GKE Autopilot] ``` +## The triage model + +A **case** is a proband: a VCF plus the HPO terms observed in that patient. Annotation produces +variants; the model scores them; ranking then answers the only question that matters — which few +variants could explain *this* phenotype. + +Rarity (<0.1% in gnomAD) and consequence (HIGH or MODERATE) *filter*, which is the usual first +pass. Phenotype only *ranks*: a real diagnosis can sit in a gene nobody has annotated yet, and +filtering on phenotype would hide exactly that case. The rank is a weighted sum whose parts are +shown next to every candidate (`app/services/triage.py`): + +| Component | Weight | +|---|---| +| phenotype terms of this patient annotated to the gene | 0.35 | +| rarity in gnomAD | 0.25 | +| consequence severity | 0.20 | +| model P(pathogenic) | 0.20 | + +**ClinVar is deliberately not an input.** It sits beside the result as independent confirmation, so +the demo never ranks a variant highly merely because ClinVar already called it pathogenic. On the +simulated NF2 case the planted variant ranks first on phenotype, rarity and consequence alone, and +ClinVar agrees afterwards. + +Each candidate can be shortlisted or dismissed with a reason and a note; the case report is that +decision trail plus the funnel counts and the provenance (VEP version, model version, run time). +There is no authentication, so decisions are shared by everyone who opens the demo. + ## Two deployment tracks The same images and the same pipeline, deployed two ways (`infra/terraform/variables.tf`): diff --git a/docs/data.md b/docs/data.md index a9d0925..be67868 100644 --- a/docs/data.md +++ b/docs/data.md @@ -28,6 +28,7 @@ benchmarks variant callers against, so it is both realistic and unambiguously sh | **gnomAD** v4 | allele frequency feature and filter | `gs://gcp-public-data--gnomad`, `s3://gnomad-public-us-east-1` | free use, no restriction | Chen et al., *Nature* 625:92–100, 2024. [10.1038/s41586-023-06045-0](https://doi.org/10.1038/s41586-023-06045-0); Karczewski et al., *Nature* 581:434–443, 2020. [10.1038/s41586-020-2308-7](https://doi.org/10.1038/s41586-020-2308-7) | | **1000 Genomes** 30x | optional cohort/trio data | EBI FTP, `s3://1000genomes` | fully open, no access restriction | Byrska-Bishop et al., *Cell* 185(18):3426–3440.e19, 2022. [10.1016/j.cell.2022.08.004](https://doi.org/10.1016/j.cell.2022.08.004) | | **MANE Select** | one transcript per gene, if transcript choice ever matters | Ensembl/RefSeq | open | Morales et al., *Nature* 604:310–315, 2022. [10.1038/s41586-022-04558-8](https://doi.org/10.1038/s41586-022-04558-8) | +| **Human Phenotype Ontology** gene-to-phenotype | what the phenotype half of the ranking matches against (`make hpo`) | `purl.obolibrary.org/obo/hp/hpoa/genes_to_phenotype.txt` | free to use with attribution | Gargano et al., *Nucleic Acids Res* 52(D1):D1333–D1346, 2024. [10.1093/nar/gkad1005](https://doi.org/10.1093/nar/gkad1005) | ## Tools and scores @@ -40,6 +41,24 @@ benchmarks variant callers against, so it is both realistic and unambiguously sh Neither score is required: `rarelens_ml.features` treats a missing CADD or AlphaMissense value as NaN and LightGBM handles it, so the pipeline runs without the plugin data. +## The simulated proband + +`make demo-case` builds `data/proband-simulated.vcf.gz`: real GIAB HG002 variants as background +plus one real ClinVar 2-star pathogenic variant in a disease gene (NF2 by default, giving +neurofibromatosis type 2). Spiking a known variant into a public genome is how phenotype-driven +triage tools are benchmarked, every input is public, and the file's header says SIMULATED. It is +not a patient, and no part of it is invented: both the background and the planted variant are real +published records. + +The point of it is that the case has a right answer, so the ranking can be checked rather than +admired. Give the case the phenotype of the planted disease and the planted variant should rank +first — on phenotype, rarity and consequence, with ClinVar agreeing only afterwards. + +**Caveat when running without a VEP cache.** `VEP_DATABASE=true` queries Ensembl's public database +instead of the 25 GB cache. It returns no gnomAD frequencies, so every variant looks absent from +gnomAD and the rarity term stops discriminating. Fine for showing the mechanics; use the cache for +anything you would quote. + ## Evaluating the model honestly The model trains on ClinVar labels and is scored on ClinVar-labelled variants, which is exactly diff --git a/scripts/load-hpo.py b/scripts/load-hpo.py new file mode 100755 index 0000000..e380aeb --- /dev/null +++ b/scripts/load-hpo.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python3 +"""Load HPO's gene-to-phenotype annotations into the gene_phenotypes table. + +This is the reference data the ranking matches a case's phenotype against. Source file: +https://purl.obolibrary.org/obo/hp/hpoa/genes_to_phenotype.txt (HPO release, ~20 MB). + +Cite the Human Phenotype Ontology when showing results; see docs/data.md. +""" +import argparse +import csv +import io +import os +import sys +import urllib.request + +from sqlalchemy import create_engine, text +from sqlalchemy.engine import make_url + +URL = "https://purl.obolibrary.org/obo/hp/hpoa/genes_to_phenotype.txt" + + +def rows(handle: io.TextIOBase) -> list[tuple[str, str, str]]: + """Unique (gene, term) pairs; the file repeats them once per associated disease.""" + seen: set[tuple[str, str]] = set() + out: list[tuple[str, str, str]] = [] + for row in csv.DictReader(handle, delimiter="\t"): + gene, hpo_id, name = row["gene_symbol"], row["hpo_id"], row["hpo_name"] + if not gene or not hpo_id or (gene, hpo_id) in seen: + continue + seen.add((gene, hpo_id)) + out.append((gene[:60], hpo_id[:20], name[:200])) + return out + + +def main() -> None: + p = argparse.ArgumentParser() + p.add_argument("--url", default=URL) + p.add_argument("--file", help="use a local copy instead of downloading") + a = p.parse_args() + + url = os.environ.get("DATABASE_URL") + if not url: + sys.exit("DATABASE_URL is not set") + + if a.file: + with open(a.file) as fh: + annotations = rows(fh) + else: + print(f"downloading {a.url}", file=sys.stderr) + with urllib.request.urlopen(a.url) as response: # noqa: S310 - fixed HPO release URL + annotations = rows(io.TextIOWrapper(response, encoding="utf-8")) + print(f"{len(annotations)} gene/term pairs", file=sys.stderr) + + engine = create_engine(make_url(url).set(drivername="postgresql+psycopg")) + with engine.begin() as conn: + conn.execute(text("TRUNCATE gene_phenotypes RESTART IDENTITY")) + cursor = conn.connection.cursor() + with cursor.copy("COPY gene_phenotypes (gene_symbol, hpo_id, hpo_name) FROM STDIN") as copy: + for row in annotations: + copy.write_row(row) + print(f"loaded {len(annotations)} annotations", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/scripts/make-demo-case.sh b/scripts/make-demo-case.sh new file mode 100755 index 0000000..ebc382e --- /dev/null +++ b/scripts/make-demo-case.sh @@ -0,0 +1,52 @@ +#!/usr/bin/env bash +# Build a SIMULATED proband VCF: real GIAB HG002 variants as background, plus one real ClinVar +# pathogenic variant in a disease gene, which becomes the diagnosis the triage should find. +# +# Spiking a known variant into a public genome is how phenotype-driven triage tools are +# benchmarked. Every input here is public and openly licensed, and this is not a real patient. +# Provenance and citations: docs/data.md +set -euo pipefail + +IMAGE=${BCFTOOLS_IMAGE:-quay.io/biocontainers/bcftools:1.20--h8b25389_0} +GENE=${GENE:-NF2} +REGION=${REGION:-chr22:20000000-31000000} +BACKGROUND=${BACKGROUND:-12} # kept small: VEP's database mode is slow per variant +OUT_DIR=${OUT_DIR:-data} +CLINVAR=${CLINVAR:-https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz} +GIAB=${GIAB:-https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NISTv4.2.1/GRCh38/HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz} + +mkdir -p "$OUT_DIR" +docker run --rm -v "$PWD/$OUT_DIR:/out" "$IMAGE" bash -eu -c " + echo '==> background: GIAB HG002 $REGION' >&2 + bcftools view -H -r '$REGION' '$GIAB' \ + | awk '\$5 !~ /,/ && length(\$4) < 20 && length(\$5) < 20' \ + | awk 'NR % 149 == 0' \ + | head -n $BACKGROUND \ + | awk -v OFS='\t' '{ sub(/^chr/, \"\", \$1); print \$1, \$2, \".\", \$4, \$5, \".\", \"PASS\", \".\" }' \ + > /out/_background.tsv + + echo '==> diagnosis: a ClinVar pathogenic variant in $GENE' >&2 + bcftools view -H -r 22 -i 'INFO/CLNSIG ~ \"Pathogenic\" && INFO/GENEINFO ~ \"${GENE}:\" && INFO/CLNREVSTAT ~ \"multiple_submitters\"' '$CLINVAR' \ + | awk 'length(\$4) < 20 && length(\$5) < 20' \ + | head -n 1 \ + | awk -v OFS='\t' '{ print \$1, \$2, \".\", \$4, \$5, \".\", \"PASS\", \".\" }' \ + > /out/_spike.tsv + + test -s /out/_background.tsv || { echo 'no background variants found' >&2; exit 1; } + test -s /out/_spike.tsv || { echo 'no 2-star pathogenic ClinVar variant found for $GENE' >&2; exit 1; } + + { + printf '##fileformat=VCFv4.2\n' + printf '##source=rarelens SIMULATED proband: GIAB HG002 background + one ClinVar pathogenic %s variant\n' '$GENE' + printf '##contig=\n' + printf '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n' + sort -k2,2n /out/_background.tsv /out/_spike.tsv + } | bgzip > /out/proband-simulated.vcf.gz + tabix -f -p vcf /out/proband-simulated.vcf.gz + rm -f /out/_background.tsv +" +echo +echo "wrote $OUT_DIR/proband-simulated.vcf.gz" +echo "the planted diagnosis (chrom pos ref alt):" +awk '{print " " $1, $2, $4, $5}' "$OUT_DIR/_spike.tsv" +rm -f "$OUT_DIR/_spike.tsv" diff --git a/web/src/app.css b/web/src/app.css index 070f438..45f4e73 100644 --- a/web/src/app.css +++ b/web/src/app.css @@ -12,7 +12,7 @@ } html { background: var(--paper); color: var(--ink); font-family: var(--sans); font-size: 17px; } body { margin: 0; } -main { max-width: 1100px; margin: 0 auto; padding: 2.5rem 1.5rem; } +main { max-width: 1280px; margin: 0 auto; padding: 2.5rem 1.5rem; } h1 { font-weight: 600; font-size: 2rem; letter-spacing: -0.01em; margin: 0 0 0.5rem; } h2 { font-weight: 600; font-size: 1.25rem; margin: 2rem 0 0.75rem; } p.lede { color: var(--ink-soft); max-width: 60ch; margin: 0 0 2rem; } @@ -40,3 +40,89 @@ td.coord, td.hgvs { font-family: var(--mono); font-size: 0.85rem; } /* aligned .spinner { animation: spin 0.9s linear infinite; } } @keyframes spin { to { transform: rotate(360deg); } } + +/* --- triage ------------------------------------------------------------------------------- */ +.muted { color: var(--ink-soft); } +.actions { display: flex; align-items: center; gap: 1rem; } +.reportlink { font-weight: 600; } +.provenance { font-size: 0.85rem; margin: 0.35rem 0 0; } +.casehead { margin-bottom: 1.25rem; } +.casehead h1 { margin-bottom: 0.4rem; } +.log { white-space: pre-wrap; background: white; border: 1px solid var(--line); padding: 0.75rem; + max-height: 14rem; overflow: auto; font-size: 0.8rem; } + +.chips { display: inline-flex; flex-wrap: wrap; gap: 0.35rem; align-items: center; } +.chip { font-size: 0.78rem; padding: 0.12rem 0.5rem; border-radius: 999px; border: 1px solid var(--line); + background: white; color: var(--ink-soft); white-space: nowrap; } +.chip.match { background: var(--plum-soft); border-color: transparent; color: var(--plum); font-weight: 600; } +.chip.rare { border-color: var(--ink-soft); color: var(--ink); } +.chip.impact { background: var(--ink); border-color: var(--ink); color: white; } +.chip.model { border-style: dashed; color: var(--ink); } +.chip.clinvar { background: #fdf1df; border-color: transparent; color: var(--amber); font-weight: 600; } +.chipx { background: none; border: none; color: inherit; padding: 0 0 0 0.3rem; cursor: pointer; font: inherit; } + +/* The funnel is the headline: thousands of variants down to a handful. */ +.funnel { list-style: none; padding: 0; margin: 0 0 1.5rem; display: grid; gap: 0.35rem; } +.funnel li { display: grid; grid-template-columns: 5rem 16rem 1fr; align-items: center; gap: 0.75rem; } +.funnel-value { font-family: var(--mono); font-size: 1.05rem; text-align: right; font-variant-numeric: tabular-nums; } +.funnel-label { color: var(--ink-soft); font-size: 0.9rem; } +.funnel-track { background: white; border: 1px solid var(--line); height: 0.75rem; border-radius: 2px; } +.funnel-bar { display: block; height: 100%; background: var(--plum); opacity: 0.75; } + +.filters { display: flex; gap: 0.5rem; flex-wrap: wrap; align-items: center; margin-bottom: 1rem; } +.triage { display: grid; grid-template-columns: minmax(0, 1fr); gap: 1.25rem; align-items: start; } +@media (min-width: 62rem) { .triage { grid-template-columns: minmax(0, 1fr) 24rem; } } + +.candidate { display: grid; grid-template-columns: 2rem minmax(10rem, 16rem) 1fr auto auto; gap: 0.75rem; + align-items: center; width: 100%; text-align: left; background: white; color: var(--ink); + border: 1px solid var(--line); border-left: 3px solid transparent; border-radius: 4px; + padding: 0.6rem 0.75rem; margin-bottom: 0.4rem; font: inherit; cursor: pointer; } +.candidate:hover { border-left-color: var(--plum-soft); } +.candidate.selected { border-left-color: var(--plum); background: #fbf7fa; } +.candidate .rank { font-family: var(--mono); color: var(--ink-soft); } +.candidate .who { display: flex; flex-direction: column; gap: 0.1rem; } +.candidate .who .coord, .candidate .who .hgvs { font-size: 0.78rem; color: var(--ink-soft); } +.rankscore { display: flex; align-items: center; gap: 0.5rem; } +.scorebar { width: 4rem; height: 0.5rem; background: var(--line); border-radius: 2px; overflow: hidden; } +.scorebar span { display: block; height: 100%; background: var(--plum); } +.scorenum { font-family: var(--mono); font-size: 0.85rem; font-variant-numeric: tabular-nums; } +.decision { font-size: 0.72rem; text-transform: uppercase; letter-spacing: 0.04em; font-weight: 600; } +.decision.shortlisted { color: var(--plum); } +.decision.dismissed { color: var(--ink-soft); } + +.panel { background: white; border: 1px solid var(--line); border-radius: 4px; padding: 1rem; position: sticky; top: 1rem; } +.panel-head { display: flex; justify-content: space-between; align-items: start; gap: 0.5rem; } +.panel h3 { margin: 0 0 0.5rem; font-size: 1.05rem; } +.panel h4 { margin: 1.1rem 0 0.4rem; font-size: 0.78rem; text-transform: uppercase; + letter-spacing: 0.05em; color: var(--ink-soft); } +.panel input, .panel textarea { width: 100%; margin-bottom: 0.5rem; font-family: inherit; } +.components { border: none; background: none; } +.components th { font-weight: 400; color: var(--ink-soft); padding: 0.15rem 0.5rem 0.15rem 0; } +.components td { padding: 0.15rem 0; font-size: 0.82rem; text-align: right; } +.evidence { display: grid; grid-template-columns: 7rem 1fr; gap: 0.2rem 0.5rem; margin: 0; font-size: 0.88rem; } +.evidence dt { color: var(--ink-soft); } +.evidence dd { margin: 0; } +.links { display: flex; gap: 0.75rem; font-size: 0.85rem; } +.decide { display: flex; gap: 0.5rem; } + +.newcase { display: grid; gap: 0.6rem; margin-bottom: 1.5rem; } +.newcase .row { display: flex; gap: 0.5rem; flex-wrap: wrap; } +.newcase button { justify-self: start; } +.picker { position: relative; display: grid; gap: 0.4rem; } +.picker .results { list-style: none; margin: 0; padding: 0.25rem; background: white; + border: 1px solid var(--line); border-radius: 4px; max-height: 12rem; overflow: auto; } +.picker .results button { width: 100%; text-align: left; border: none; padding: 0.3rem 0.4rem; } +.picker .results button:hover { background: var(--plum-soft); } + +.caselist { list-style: none; padding: 0; margin: 0; } +.caselist li { display: flex; flex-wrap: wrap; gap: 0.75rem; align-items: center; padding: 0.7rem 0.75rem; + background: white; border: 1px solid var(--line); border-radius: 4px; margin-bottom: 0.4rem; } +.casename { font-weight: 600; min-width: 10rem; } + +.report { background: white; border: 1px solid var(--line); padding: 2rem; } +.reportitem { border-top: 1px solid var(--line); padding-top: 0.75rem; margin-top: 0.75rem; } +.reportitem h3 { margin: 0 0 0.4rem; } +.dismissed { color: var(--ink-soft); } +.disclaimer { margin-top: 2rem; padding-top: 1rem; border-top: 1px solid var(--line); + color: var(--ink-soft); font-size: 0.85rem; } +@media print { .noprint, nav { display: none; } .report { border: none; padding: 0; } } diff --git a/web/src/lib/api.ts b/web/src/lib/api.ts index 15fbcaa..2009825 100644 --- a/web/src/lib/api.ts +++ b/web/src/lib/api.ts @@ -1,22 +1,71 @@ // Dynamic, not static: the same image serves /api behind the ingress and a full URL elsewhere. import { env } from '$env/dynamic/public'; -import { variantQuery, type VariantFilters } from './query'; +import { candidateQuery, type CandidateFilters } from './candidates'; export type Assembly = 'GRCh38' | 'GRCh37'; -export type Sample = { id: string; name: string; vcf_uri: string; assembly: Assembly; created_at: string }; +export type JobStatus = 'queued' | 'running' | 'succeeded' | 'failed'; +export type PhenotypeTerm = { hpo_id: string; label: string }; + export type Job = { - id: string; sample_id: string; status: 'queued' | 'running' | 'succeeded' | 'failed'; - log: string | null; created_at: string; finished_at: string | null; + id: string; case_id: string; status: JobStatus; + vep_version: string | null; log: string | null; created_at: string; finished_at: string | null; }; + +export type Case = { + id: string; name: string; vcf_uri: string; assembly: Assembly; created_at: string; + phenotypes: PhenotypeTerm[]; latest_job: Job | null; shortlisted: number; +}; + export type Variant = { id: number; chrom: string; pos: number; ref: string; alt: string; gene: string | null; consequence: string | null; impact: string | null; hgvsc: string | null; hgvsp: string | null; - gnomad_af: number | null; clinvar_sig: string | null; prediction: { score: number; model_version: string } | null; + gnomad_af: number | null; clinvar_sig: string | null; + prediction: { score: number; model_name: string; model_version: string } | null; +}; + +export type DecisionState = 'shortlisted' | 'dismissed'; +export type Decision = { + state: DecisionState; reason: string | null; note: string | null; decided_at: string; +}; + +export type Candidate = { + variant: Variant; + score: number; + components: Record; + matched_terms: PhenotypeTerm[]; + scored: boolean; + decision: Decision | null; +}; + +export type VariantDetail = Candidate & { annotations: Record }; +export type Funnel = { total: number; rare: number; candidates: number; phenotype_matched: number }; + +export type CandidatePage = { + funnel: Funnel; + weights: Record; + items: Candidate[]; + total: number; + limit: number; + offset: number; +}; + +export type Report = { + case: Case; + funnel: Funnel; + generated_at: string; + provenance: { + job_id: string | null; vep_version: string | null; finished_at: string | null; + model_name: string | null; model_version: string | null; + }; + shortlisted: Candidate[]; + dismissed: Candidate[]; }; -export type VariantPage = { items: Variant[]; total: number; limit: number; offset: number }; async function req(path: string, init?: RequestInit): Promise { - const r = await fetch(`${env.PUBLIC_API_URL ?? '/api'}${path}`, { headers: { 'content-type': 'application/json' }, ...init }); + const r = await fetch(`${env.PUBLIC_API_URL ?? '/api'}${path}`, { + headers: { 'content-type': 'application/json' }, + ...init + }); if (!r.ok) throw new Error(await errorReason(r)); return r.json() as Promise; } @@ -34,12 +83,26 @@ async function errorReason(r: Response): Promise { return fallback; } +const post = (body?: unknown): RequestInit => ({ + method: 'POST', + body: body === undefined ? undefined : JSON.stringify(body) +}); + export const api = { - samples: () => req('/samples'), - createSample: (body: Pick) => - req('/samples', { method: 'POST', body: JSON.stringify(body) }), - jobsForSample: (sampleId: string) => req(`/samples/${sampleId}/jobs`), - annotate: (sampleId: string) => req(`/samples/${sampleId}/annotate`, { method: 'POST' }), + cases: () => req('/cases'), + createCase: (body: { name: string; vcf_uri: string; assembly: Assembly; phenotypes: PhenotypeTerm[] }) => + req('/cases', post(body)), + case: (id: string) => req(`/cases/${id}`), + annotate: (id: string) => req(`/cases/${id}/annotate`, post()), + score: (id: string) => req<{ scored: number; model_version: string }>(`/cases/${id}/score`, post()), job: (id: string) => req(`/jobs/${id}`), - variants: (filters: VariantFilters) => req(`/variants?${variantQuery(filters)}`) + candidates: (caseId: string, filters: CandidateFilters = {}) => { + const query = candidateQuery(filters); + return req(`/cases/${caseId}/candidates${query ? `?${query}` : ''}`); + }, + variant: (id: number) => req(`/variants/${id}`), + decide: (id: number, body: { state: DecisionState; reason?: string; note?: string }) => + req(`/variants/${id}/decision`, post(body)), + report: (caseId: string) => req(`/cases/${caseId}/report`), + phenotypes: (q: string) => req(`/phenotypes?q=${encodeURIComponent(q)}`) }; diff --git a/web/src/lib/candidates.test.ts b/web/src/lib/candidates.test.ts new file mode 100644 index 0000000..e75e664 --- /dev/null +++ b/web/src/lib/candidates.test.ts @@ -0,0 +1,88 @@ +import { describe, expect, it } from 'vitest'; +import { candidateQuery, evidenceChips, funnelSteps, scoreBarPercent } from './candidates'; +import type { Candidate, Funnel } from './api'; + +const candidate = (over: Partial = {}): Candidate => ({ + variant: { + id: 1, chrom: '22', pos: 100, ref: 'A', alt: 'G', gene: 'NF2', + consequence: 'frameshift_variant', impact: 'HIGH', hgvsc: null, hgvsp: 'p.Ser486fs', + gnomad_af: null, clinvar_sig: null, prediction: null + }, + score: 0.88, + components: { phenotype: 0.67, rarity: 1, consequence: 1, model: 0.94 }, + matched_terms: [{ hpo_id: 'HP:0000365', label: 'Hearing impairment' }], + scored: true, + decision: null, + ...over +}); + +describe('evidenceChips', () => { + it('leads with the phenotype match, because that is what makes it a candidate', () => { + const [first] = evidenceChips(candidate(), 3); + expect(first).toEqual({ label: '1/3 phenotype terms', tone: 'match' }); + }); + + it('says plainly when the gene has no link to the phenotype', () => { + const chips = evidenceChips(candidate({ matched_terms: [] }), 3); + expect(chips[0]).toEqual({ label: 'no phenotype match', tone: 'muted' }); + }); + + it('calls an absent variant absent rather than showing a zero', () => { + expect(evidenceChips(candidate(), 3)[1]).toEqual({ label: 'absent from gnomAD', tone: 'rare' }); + }); + + it('shows a frequency when there is one', () => { + const c = candidate(); + c.variant.gnomad_af = 0.00042; + expect(evidenceChips(c, 3)[1].label).toBe('gnomAD 4.2e-4'); + }); + + it('reads the consequence without the underscores', () => { + expect(evidenceChips(candidate(), 3)[2]).toEqual({ label: 'frameshift variant', tone: 'impact' }); + }); + + it('marks an unscored variant instead of implying a zero score', () => { + const chips = evidenceChips(candidate({ scored: false }), 3); + expect(chips[3]).toEqual({ label: 'unscored', tone: 'muted' }); + }); + + it('shows ClinVar last, as confirmation rather than an input', () => { + const c = candidate(); + c.variant.clinvar_sig = 'pathogenic'; + expect(evidenceChips(c, 3).at(-1)).toEqual({ label: 'ClinVar: pathogenic', tone: 'clinvar' }); + }); +}); + +describe('funnelSteps', () => { + it('describes each narrowing step in order', () => { + const funnel: Funnel = { total: 1284, rare: 41, candidates: 12, phenotype_matched: 6 }; + expect(funnelSteps(funnel).map((s) => `${s.label} ${s.value}`)).toEqual([ + 'variants called 1284', + 'rare (<0.1%) 41', + 'coding candidates 12', + 'in phenotype-matched genes 6' + ]); + }); +}); + +describe('scoreBarPercent', () => { + it('maps the 0..1 score onto a bar width', () => { + expect(scoreBarPercent(0.884)).toBe(88); + expect(scoreBarPercent(0)).toBe(0); + expect(scoreBarPercent(1)).toBe(100); + }); +}); + +describe('candidateQuery', () => { + it('sends only the filters that are set', () => { + expect(candidateQuery({})).toBe(''); + expect(candidateQuery({ gene: ' nf2 ', impact: 'HIGH', state: 'undecided' })).toBe( + 'gene=nf2&impact=HIGH&state=undecided' + ); + }); + + it('drops a cleared allele frequency box', () => { + expect(candidateQuery({ maxAf: null })).toBe(''); + expect(candidateQuery({ maxAf: 0.01 })).toBe('max_af=0.01'); + }); +}); diff --git a/web/src/lib/candidates.ts b/web/src/lib/candidates.ts new file mode 100644 index 0000000..f395c6c --- /dev/null +++ b/web/src/lib/candidates.ts @@ -0,0 +1,70 @@ +import type { Candidate, Funnel } from './api'; + +export type Chip = { + label: string; + tone: 'match' | 'rare' | 'impact' | 'model' | 'clinvar' | 'muted'; +}; + +export type CandidateFilters = { + gene?: string; + impact?: string; + maxAf?: number | string | null; + state?: 'shortlisted' | 'dismissed' | 'undecided' | ''; + limit?: number; + offset?: number; +}; + +const formatAf = (af: number) => af.toExponential(1); + +/** + * The reasons this variant is a candidate, in the order a reviewer reads them. ClinVar comes last + * and is styled apart, because it confirms the ranking rather than feeding it. + */ +export function evidenceChips(candidate: Candidate, caseTermCount: number): Chip[] { + const { variant } = candidate; + const chips: Chip[] = [ + candidate.matched_terms.length + ? { label: `${candidate.matched_terms.length}/${caseTermCount} phenotype terms`, tone: 'match' } + : { label: 'no phenotype match', tone: 'muted' }, + variant.gnomad_af === null + ? { label: 'absent from gnomAD', tone: 'rare' } + : { label: `gnomAD ${formatAf(variant.gnomad_af)}`, tone: 'rare' }, + { + label: (variant.consequence ?? 'unknown consequence').split('&')[0].replace(/_/g, ' '), + tone: 'impact' + }, + candidate.scored && variant.prediction + ? { label: `model ${variant.prediction.score.toFixed(2)}`, tone: 'model' } + : { label: 'unscored', tone: 'muted' } + ]; + if (variant.clinvar_sig) chips.push({ label: `ClinVar: ${variant.clinvar_sig}`, tone: 'clinvar' }); + return chips; +} + +/** Thousands of variants down to a handful: the narrowing is the story. */ +export function funnelSteps(funnel: Funnel): { label: string; value: number }[] { + return [ + { label: 'variants called', value: funnel.total }, + { label: 'rare (<0.1%)', value: funnel.rare }, + { label: 'coding candidates', value: funnel.candidates }, + { label: 'in phenotype-matched genes', value: funnel.phenotype_matched } + ]; +} + +export const scoreBarPercent = (score: number): number => + Math.max(0, Math.min(100, Math.round(score * 100))); + +export function candidateQuery(filters: CandidateFilters): string { + const q = new URLSearchParams(); + const gene = filters.gene?.trim(); + if (gene) q.set('gene', gene); + if (filters.impact) q.set('impact', filters.impact); + const { maxAf } = filters; + if (maxAf !== null && maxAf !== undefined && maxAf !== '' && Number.isFinite(Number(maxAf))) { + q.set('max_af', String(maxAf)); + } + if (filters.state) q.set('state', filters.state); + if (filters.limit !== undefined) q.set('limit', String(filters.limit)); + if (filters.offset !== undefined) q.set('offset', String(filters.offset)); + return q.toString(); +} diff --git a/web/src/lib/components/CandidateRow.svelte b/web/src/lib/components/CandidateRow.svelte new file mode 100644 index 0000000..121a41f --- /dev/null +++ b/web/src/lib/components/CandidateRow.svelte @@ -0,0 +1,38 @@ + + + diff --git a/web/src/lib/components/Chips.svelte b/web/src/lib/components/Chips.svelte new file mode 100644 index 0000000..e4a1810 --- /dev/null +++ b/web/src/lib/components/Chips.svelte @@ -0,0 +1,10 @@ + + + + {#each chips as chip (chip.label)} + {chip.label} + {/each} + diff --git a/web/src/lib/components/Funnel.svelte b/web/src/lib/components/Funnel.svelte new file mode 100644 index 0000000..d601ddc --- /dev/null +++ b/web/src/lib/components/Funnel.svelte @@ -0,0 +1,18 @@ + + +
    + {#each steps as step (step.label)} +
  1. + {step.value.toLocaleString('en-GB')} + {step.label} + +
    2. + {/each} +
diff --git a/web/src/lib/components/PhenotypePicker.svelte b/web/src/lib/components/PhenotypePicker.svelte new file mode 100644 index 0000000..46f9f73 --- /dev/null +++ b/web/src/lib/components/PhenotypePicker.svelte @@ -0,0 +1,62 @@ + + +
+ + {#if searching}searching…{/if} + {#if results.length} +
    + {#each results as term (term.hpo_id)} +
  • + +
    • + {/each} +
+ {/if} + + {#each selected as term (term.hpo_id)} + + {term.label} + + + {/each} + +
diff --git a/web/src/lib/components/VariantPanel.svelte b/web/src/lib/components/VariantPanel.svelte new file mode 100644 index 0000000..69059cc --- /dev/null +++ b/web/src/lib/components/VariantPanel.svelte @@ -0,0 +1,111 @@ + + + diff --git a/web/src/lib/query.test.ts b/web/src/lib/query.test.ts deleted file mode 100644 index dd602d2..0000000 --- a/web/src/lib/query.test.ts +++ /dev/null @@ -1,29 +0,0 @@ -import { describe, expect, it } from 'vitest'; -import { variantQuery } from './query'; - -describe('variantQuery', () => { - it('always includes the job id', () => { - expect(variantQuery({ jobId: 'abc' })).toBe('job_id=abc'); - }); - - it('includes filters that have values', () => { - const q = new URLSearchParams( - variantQuery({ jobId: 'abc', gene: ' brca1 ', impact: 'HIGH', maxAf: 0.01, limit: 100 }) - ); - expect(q.get('gene')).toBe('brca1'); - expect(q.get('impact')).toBe('HIGH'); - expect(q.get('max_af')).toBe('0.01'); - expect(q.get('limit')).toBe('100'); - }); - - it('drops a cleared number input instead of sending "null"', () => { - // Svelte binds an emptied to null. - for (const maxAf of [null, undefined, '', Number.NaN]) { - expect(variantQuery({ jobId: 'abc', maxAf })).toBe('job_id=abc'); - } - }); - - it('drops blank gene and impact', () => { - expect(variantQuery({ jobId: 'abc', gene: ' ', impact: '' })).toBe('job_id=abc'); - }); -}); diff --git a/web/src/lib/query.ts b/web/src/lib/query.ts deleted file mode 100644 index 1d83f74..0000000 --- a/web/src/lib/query.ts +++ /dev/null @@ -1,22 +0,0 @@ -export type VariantFilters = { - jobId: string; - gene?: string; - impact?: string; - // An emptied binds to null, so accept it and drop it. - maxAf?: number | string | null; - limit?: number; - offset?: number; -}; - -export function variantQuery(f: VariantFilters): string { - const q = new URLSearchParams({ job_id: f.jobId }); - const gene = f.gene?.trim(); - if (gene) q.set('gene', gene); - if (f.impact) q.set('impact', f.impact); - if (f.maxAf !== null && f.maxAf !== undefined && f.maxAf !== '' && Number.isFinite(Number(f.maxAf))) { - q.set('max_af', String(f.maxAf)); - } - if (f.limit !== undefined) q.set('limit', String(f.limit)); - if (f.offset !== undefined) q.set('offset', String(f.offset)); - return q.toString(); -} diff --git a/web/src/routes/+page.svelte b/web/src/routes/+page.svelte index 7b258a0..4286529 100644 --- a/web/src/routes/+page.svelte +++ b/web/src/routes/+page.svelte @@ -1,18 +1,20 @@

rarelens

-

Annotate a VCF with Ensembl VEP, score each variant, and browse what came back. Public test data only.

+

+ Rare disease triage on public data: a proband's variants narrowed against their phenotype, each + candidate carrying the evidence for its rank. Research demo, not a diagnostic tool. +

-

Add a sample

-
- - - - +

New case

+
+
+ + + +
+ +
-{#if error}

Could not add the sample: {error}

{/if} +{#if error}

Could not create the case: {error}

{/if} -

Samples

+

Cases

{#if loadError} -

Could not load samples: {loadError}.

+

Could not load cases: {loadError}.

{:else if !loaded} -

Loading samples…

-{:else if samples.length === 0} -
No samples yet. Add one above to run the annotation pipeline.
+

Loading…

+{:else if cases.length === 0} +
No cases yet. Create one above, then run the annotation pipeline on it.
{:else} - - - - {#each samples as s (s.id)} - - - - - - - {/each} - -
NameVCFAssemblyAdded
{s.name}{s.vcf_uri}{s.assembly}{new Date(s.created_at).toLocaleDateString('en-GB')}
+
    + {#each cases as c (c.id)} +
  • + {c.name} + + {#each c.phenotypes.slice(0, 4) as term (term.hpo_id)} + {term.label} + {/each} + {#if c.phenotypes.length > 4}+{c.phenotypes.length - 4}{/if} + {#if c.phenotypes.length === 0}no phenotype{/if} + + {status(c)} + {c.shortlisted} shortlisted +
    • + {/each} +
{/if} diff --git a/web/src/routes/cases/[id]/+page.svelte b/web/src/routes/cases/[id]/+page.svelte new file mode 100644 index 0000000..5be56b3 --- /dev/null +++ b/web/src/routes/cases/[id]/+page.svelte @@ -0,0 +1,206 @@ + + +All cases +
+

{kase?.name ?? 'Case'}

+ + {#each kase?.phenotypes ?? [] as term (term.hpo_id)} + {term.label} + {/each} + {#if kase && kase.phenotypes.length === 0} + no phenotype recorded — ranking falls back to rarity and consequence + {/if} + +

+ {kase?.assembly} + {#if job?.vep_version}· VEP {job.vep_version}{/if} + {#if page?.items[0]?.variant.prediction}· model {page.items[0].variant.prediction.model_version}{/if} + {#if job?.finished_at}· analysed {new Date(job.finished_at).toLocaleString('en-GB')}{/if} +

+
+ +{#if error}

{error}

{/if} + +{#if running} +

+ + Analysing for {formatElapsed(elapsedMs)}{#if step} · {step}{/if} +

+{:else} +

+ + {#if page && page.funnel.total > 0} + Case report → + {/if} +

+{/if} + +{#if job?.status === 'failed' && job.log} + 
{job.log}
+{/if} + +{#if page && page.funnel.total > 0} + + +
+ + + + + +
+ +
+
+

{page.total} candidates ranked by phenotype fit, rarity, consequence and model score

+ {#each page.items as candidate, i (candidate.variant.id)} + select(candidate)} /> + {:else} +
No candidates match these filters.
+ {/each} +
+ {#if detail && page} + (detail = null)} + ondecided={decided} /> + {/if} +
+{:else if page && !running} +
+ No results yet. "Analyse case" runs VEP over the VCF, scores each variant, then ranks what is + left against the phenotype. +
+{/if} diff --git a/web/src/routes/cases/[id]/+page.ts b/web/src/routes/cases/[id]/+page.ts new file mode 100644 index 0000000..2b4618b --- /dev/null +++ b/web/src/routes/cases/[id]/+page.ts @@ -0,0 +1,2 @@ +import type { PageLoad } from './$types'; +export const load: PageLoad = ({ params }) => ({ caseId: params.id }); diff --git a/web/src/routes/cases/[id]/report/+page.svelte b/web/src/routes/cases/[id]/report/+page.svelte new file mode 100644 index 0000000..3a5211c --- /dev/null +++ b/web/src/routes/cases/[id]/report/+page.svelte @@ -0,0 +1,78 @@ + + +← Back to triage +{#if error}

{error}

{/if} + +{#if report} +
+

Case report — {report.case.name}

+

+ Generated {new Date(report.generated_at).toLocaleString('en-GB')} · + {report.case.assembly} + {#if report.provenance.vep_version}· VEP {report.provenance.vep_version}{/if} + {#if report.provenance.model_version}· model {report.provenance.model_version}{/if} +

+ +

Phenotype

+ + {#each report.case.phenotypes as term (term.hpo_id)} + {term.label} {term.hpo_id} + {:else} + none recorded + {/each} + + +

Narrowing

+ + +

Shortlisted ({report.shortlisted.length})

+ {#each report.shortlisted as item (item.variant.id)} +
+

{item.variant.gene ?? 'intergenic'} {item.variant.chrom}:{item.variant.pos} {item.variant.ref}>{item.variant.alt}

+ +

Reason: {item.decision?.reason ?? '—'}

+ {#if item.decision?.note}

{item.decision.note}

{/if} +
+ {:else} +

Nothing shortlisted yet.

+ {/each} + + {#if report.dismissed.length} +

Dismissed ({report.dismissed.length})

+
    + {#each report.dismissed as item (item.variant.id)} +
  • + {item.variant.gene ?? 'intergenic'} + {item.variant.chrom}:{item.variant.pos} + — {item.decision?.reason ?? 'no reason given'} +
    • + {/each} +
+ {/if} + +

+ Research demonstration on public data. Computational evidence is supporting only under + ACMG/AMP guidance and this tool makes no diagnostic claim. +

+
+{/if} diff --git a/web/src/routes/cases/[id]/report/+page.ts b/web/src/routes/cases/[id]/report/+page.ts new file mode 100644 index 0000000..2b4618b --- /dev/null +++ b/web/src/routes/cases/[id]/report/+page.ts @@ -0,0 +1,2 @@ +import type { PageLoad } from './$types'; +export const load: PageLoad = ({ params }) => ({ caseId: params.id }); diff --git a/web/src/routes/samples/[id]/+page.svelte b/web/src/routes/samples/[id]/+page.svelte deleted file mode 100644 index c3e3b9c..0000000 --- a/web/src/routes/samples/[id]/+page.svelte +++ /dev/null @@ -1,141 +0,0 @@ - - -All samples -

Sample {data.sampleId.slice(0, 8)}

- -{#if error}

{error}

{/if} - -{#if job} -

Job {job.id.slice(0, 8)}: {job.status}

- {#if running} -

- - Annotating for {formatElapsed(elapsedMs)}{#if step} · {step}{/if} -

- {/if} - {#if job.status === 'failed' && job.log} - 
{job.log}
- {/if} -{/if} -{#if !job || job.status === 'failed'} - -{/if} - -{#if job?.status === 'succeeded'} -

Variants

-
- - - - -
- - {#if page && page.items.length === 0} -
No variants match these filters. Raise the allele frequency cap or clear the gene filter.
- {:else if page} -

{page.total} variants, showing {page.items.length}

- - - - {#each page.items as v (v.id)} - - - - - - - - - - {/each} - -
PositionGeneConsequenceHGVSgnomAD AFClinVarScore
{v.chrom}:{v.pos} {v.ref}>{v.alt}{v.gene ?? ''}{v.consequence ?? ''}
{v.impact ?? ''}		{v.hgvsp ?? v.hgvsc ?? ''}{v.gnomad_af?.toExponential(2) ?? 'absent'}{v.clinvar_sig ?? ''}{#if v.prediction}{v.prediction.score.toFixed(2)}{:else}unscored{/if}
- {/if} -{/if} diff --git a/web/src/routes/samples/[id]/+page.ts b/web/src/routes/samples/[id]/+page.ts deleted file mode 100644 index 6effbd8..0000000 --- a/web/src/routes/samples/[id]/+page.ts +++ /dev/null @@ -1,2 +0,0 @@ -import type { PageLoad } from './$types'; -export const load: PageLoad = ({ params }) => ({ sampleId: params.id });