Samad14 commited on
Commit
7c22d1d
·
verified ·
1 Parent(s): 3cbaf1f

Upload folder using huggingface_hub

Browse files
.env.example CHANGED
@@ -21,3 +21,12 @@ DEMO_MODE=false
21
 
22
  # Hugging Face CLI (optional — for `hf upload` deploys only, not read by the app)
23
  HF_TOKEN=
 
 
 
 
 
 
 
 
 
 
21
 
22
  # Hugging Face CLI (optional — for `hf upload` deploys only, not read by the app)
23
  HF_TOKEN=
24
+
25
+ # Real nf-core/sarek executors (all fail closed when disabled or incomplete)
26
+ NGS_LOCAL_EXECUTION_ENABLED=false
27
+ NGS_SLURM_EXECUTION_ENABLED=false
28
+ NGS_AWS_BATCH_EXECUTION_ENABLED=false
29
+ NGS_RUN_ROOT=/tmp/bionexus-ngs-runs
30
+ NGS_AWS_REGION=
31
+ NGS_AWS_BATCH_JOB_QUEUE=
32
+ NGS_AWS_BATCH_JOB_DEFINITION=
NGS_PRODUCTION_EXECUTION.md ADDED
@@ -0,0 +1,50 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ # BioNexus production NGS execution
2
+
3
+ BioNexus submits the pinned `nf-core/sarek` 3.10.0 command produced by the production validator. Production execution never falls back to the exploratory Python pipeline.
4
+
5
+ ## Shared requirements
6
+
7
+ - Apply `migrations/011_ngs_production_runs.sql` to Supabase.
8
+ - Configure authentication and the existing Supabase service credentials.
9
+ - Stage the Sarek samplesheet, reference resources, optional WES BED, work directory and result directory where the selected executor can read them.
10
+ - Provide a reviewed Nextflow configuration for SLURM and AWS Batch.
11
+ - Keep the executor disabled until its preflight requirements are satisfied.
12
+
13
+ ## Local Docker, Apptainer or Singularity
14
+
15
+ Run the API on the workstation that owns the staged paths. Install Java, Nextflow, and the selected container runtime, then set:
16
+
17
+ ```env
18
+ NGS_LOCAL_EXECUTION_ENABLED=true
19
+ NGS_RUN_ROOT=/durable/bionexus-ngs-runs
20
+ ```
21
+
22
+ The detached worker writes its real Nextflow log and status under `NGS_RUN_ROOT/<run-id>/`. The output artifacts remain in the requested Sarek `--outdir`.
23
+
24
+ ## SLURM
25
+
26
+ Run the API or a private execution service on a SLURM login/head node with `nextflow`, `java`, `sbatch`, and `sacct` on `PATH`. Samples, references, work and results must use shared cluster paths.
27
+
28
+ ```env
29
+ NGS_SLURM_EXECUTION_ENABLED=true
30
+ NGS_RUN_ROOT=/shared/bionexus/ngs-runs
31
+ ```
32
+
33
+ Select the SLURM profile and provide the reviewed site-specific Nextflow configuration. BioNexus submits with `sbatch` and reads scheduler state with `sacct`.
34
+
35
+ ## AWS Batch
36
+
37
+ Create an AWS Batch queue and a driver job definition whose image includes Java, Nextflow and the AWS permissions required to submit Sarek process jobs. Grant the BioNexus API role `batch:SubmitJob`, `batch:DescribeJobs`, S3 access to the staged inputs/results, and the appropriate job-role permissions.
38
+
39
+ ```env
40
+ NGS_AWS_BATCH_EXECUTION_ENABLED=true
41
+ NGS_AWS_REGION=ap-south-1
42
+ NGS_AWS_BATCH_JOB_QUEUE=<queue-name-or-arn>
43
+ NGS_AWS_BATCH_JOB_DEFINITION=<driver-job-definition-name-or-arn>
44
+ ```
45
+
46
+ Use S3 input/result locations or mounts explicitly configured in the job definition. The artifact importer inventories the real S3 output prefix after the Batch job succeeds.
47
+
48
+ ## Scientific boundary
49
+
50
+ Successful workflow execution is not clinical validation. BioNexus imports observed execution, MultiQC, alignment, variant, coverage, identity/contamination and provenance artifacts. Missing groups remain missing and keep the clinical gate closed.
app/config.py CHANGED
@@ -32,6 +32,14 @@ class Settings:
32
  DEMO_MODE: bool = os.getenv("DEMO_MODE", "false").lower() in ("true", "1", "yes")
33
  CORS_ORIGIN: str = os.getenv("CORS_ORIGIN", "https://bioai-platform.vercel.app")
34
  SENTRY_DSN: str = os.getenv("SENTRY_DSN", "")
 
 
 
 
 
 
 
 
35
 
36
 
37
  settings = Settings()
 
32
  DEMO_MODE: bool = os.getenv("DEMO_MODE", "false").lower() in ("true", "1", "yes")
33
  CORS_ORIGIN: str = os.getenv("CORS_ORIGIN", "https://bioai-platform.vercel.app")
34
  SENTRY_DSN: str = os.getenv("SENTRY_DSN", "")
35
+ NGS_CLINICAL_EVIDENCE_HMAC_KEY: str = os.getenv("NGS_CLINICAL_EVIDENCE_HMAC_KEY", "")
36
+ NGS_LOCAL_EXECUTION_ENABLED: bool = os.getenv("NGS_LOCAL_EXECUTION_ENABLED", "false").lower() in ("true", "1", "yes")
37
+ NGS_SLURM_EXECUTION_ENABLED: bool = os.getenv("NGS_SLURM_EXECUTION_ENABLED", "false").lower() in ("true", "1", "yes")
38
+ NGS_AWS_BATCH_EXECUTION_ENABLED: bool = os.getenv("NGS_AWS_BATCH_EXECUTION_ENABLED", "false").lower() in ("true", "1", "yes")
39
+ NGS_RUN_ROOT: str = os.getenv("NGS_RUN_ROOT", "/tmp/bionexus-ngs-runs")
40
+ NGS_AWS_REGION: str = os.getenv("NGS_AWS_REGION", "")
41
+ NGS_AWS_BATCH_JOB_QUEUE: str = os.getenv("NGS_AWS_BATCH_JOB_QUEUE", "")
42
+ NGS_AWS_BATCH_JOB_DEFINITION: str = os.getenv("NGS_AWS_BATCH_JOB_DEFINITION", "")
43
 
44
 
45
  settings = Settings()
app/models/responses.py CHANGED
@@ -1,5 +1,5 @@
1
- from pydantic import BaseModel
2
- from typing import Any, Optional
3
 
4
 
5
  class PipelineRunResponse(BaseModel):
@@ -38,3 +38,96 @@ class ProfileUpdateResponse(BaseModel):
38
 
39
  class ErrorResponse(BaseModel):
40
  detail: str
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from pydantic import BaseModel, Field
2
+ from typing import Any, Literal, Optional
3
 
4
 
5
  class PipelineRunResponse(BaseModel):
 
38
 
39
  class ErrorResponse(BaseModel):
40
  detail: str
41
+
42
+
43
+ class NgsProductionPlanRequest(BaseModel):
44
+ assay: Literal["WGS", "WES"]
45
+ sample_model: Literal["singleton", "cohort", "duo", "trio", "family"] = "singleton"
46
+ input_type: Literal["FASTQ", "BAM", "CRAM"] = "FASTQ"
47
+ start_step: Literal["mapping", "markduplicates", "variant_calling"] = "mapping"
48
+ samplesheet_path: str = Field(min_length=1)
49
+ outdir: str = Field(min_length=1)
50
+ genome: Literal["GRCh38", "GRCh37"] = "GRCh38"
51
+ execution_profile: Literal["docker", "singularity", "apptainer", "slurm", "awsbatch"] = "docker"
52
+ caller: Literal["haplotypecaller", "deepvariant"] = "haplotypecaller"
53
+ target_bed: Optional[str] = None
54
+ custom_config: Optional[str] = None
55
+ annotate_with_vep: bool = True
56
+ clinical_intent: bool = False
57
+
58
+
59
+ class NgsProductionPlanResponse(BaseModel):
60
+ schema_version: str
61
+ workflow: dict[str, str]
62
+ state: Literal["PLANNED", "BLOCKED"]
63
+ ready_to_launch: bool
64
+ blockers: list[str]
65
+ warnings: list[str]
66
+ command_argv: list[str]
67
+ command_display: str
68
+ required_artifacts: list[dict[str, Any]]
69
+ provenance_requirements: list[str]
70
+ clinical_boundary: dict[str, Any]
71
+
72
+
73
+ class NgsProductionSubmitResponse(BaseModel):
74
+ run_id: str
75
+ state: Literal["SUBMITTED", "BLOCKED"]
76
+ executor: Literal["local", "slurm", "awsbatch"]
77
+ executor_job_id: Optional[str] = None
78
+ message: str
79
+
80
+
81
+ class NgsProductionRunResponse(BaseModel):
82
+ run_id: str
83
+ state: Literal["SUBMITTED", "PENDING", "RUNNING", "SUCCEEDED", "FAILED", "UNKNOWN"]
84
+ executor: Literal["local", "slurm", "awsbatch"]
85
+ executor_job_id: str
86
+ workflow: str
87
+ revision: str
88
+ outdir: str
89
+ submitted_at: str
90
+ updated_at: str
91
+ exit_code: Optional[int] = None
92
+ message: Optional[str] = None
93
+
94
+
95
+ class NgsClinicalEvidenceRequest(BaseModel):
96
+ evidence_bundle_sha256: Optional[str] = None
97
+ evidence_signature: Optional[str] = None
98
+ assay_validation_id: Optional[str] = None
99
+ workflow_status: Literal["PLANNED", "COMPLETED", "FAILED"] = "PLANNED"
100
+ sarek_revision: Optional[str] = None
101
+ reference_build: Optional[Literal["GRCh38", "GRCh37"]] = None
102
+ reference_manifest_sha256: Optional[str] = None
103
+ samplesheet_sha256: Optional[str] = None
104
+ container_digests_complete: bool = False
105
+ complete_input_processed: bool = False
106
+ required_artifacts_present: bool = False
107
+ qc_pass: bool = False
108
+ sample_identity_pass: bool = False
109
+ contamination_pass: bool = False
110
+ sex_ploidy_reviewed: bool = False
111
+ truthset_name: Optional[str] = None
112
+ benchmark_protocol_id: Optional[str] = None
113
+ benchmark_acceptance_pass: bool = False
114
+ confident_regions_sha256: Optional[str] = None
115
+ same_sample_reference_regions: bool = False
116
+ snv_precision: Optional[float] = Field(default=None, ge=0.0, le=1.0)
117
+ snv_recall: Optional[float] = Field(default=None, ge=0.0, le=1.0)
118
+ indel_precision: Optional[float] = Field(default=None, ge=0.0, le=1.0)
119
+ indel_recall: Optional[float] = Field(default=None, ge=0.0, le=1.0)
120
+ human_reviewed: bool = False
121
+ reviewer_id: Optional[str] = None
122
+ release_signature_id: Optional[str] = None
123
+ unresolved_deviations: list[str] = Field(default_factory=list)
124
+
125
+
126
+ class NgsClinicalEvidenceResponse(BaseModel):
127
+ schema_version: str
128
+ status: Literal["SOFTWARE_GATE_PASSED", "NOT_CLINICALLY_RELEASABLE"]
129
+ clinically_validated: Literal[False]
130
+ gates: list[dict[str, Any]]
131
+ missing_or_failed: list[str]
132
+ benchmark_summary: dict[str, Any]
133
+ disclaimer: str
app/ngs/contracts.py CHANGED
@@ -86,6 +86,7 @@ class ThresholdRule:
86
  evaluate: Callable[[Any], QcStatus]
87
  expectation: Optional[str] = None
88
  optional: bool = False # missing metric -> WARN (not FAIL); used for conditional sub-checks
 
89
 
90
  def apply(self, value: Any, detail: Optional[str] = None) -> Metric:
91
  return Metric(
@@ -153,7 +154,7 @@ def apply_rules(rules: list[ThresholdRule], values: dict[str, Any]) -> list[Metr
153
  elif rule.optional:
154
  metrics.append(
155
  Metric(name=rule.metric, value=None, status=QcStatus.WARN, expected=rule.expectation,
156
- detail="optional metric not evaluated (no reference data supplied)")
157
  )
158
  else:
159
  metrics.append(
@@ -227,6 +228,7 @@ class StageResult:
227
  qc: Optional[QcResult] = None
228
  decision: Decision = Decision.CONTINUE
229
  data: dict = field(default_factory=dict) # stage-specific payload (parsed results)
 
230
 
231
  def to_dict(self) -> dict:
232
  return {
@@ -241,6 +243,7 @@ class StageResult:
241
  "qc": self.qc.to_dict() if self.qc else None,
242
  "decision": self.decision.value,
243
  "data": self.data,
 
244
  }
245
 
246
 
@@ -269,6 +272,7 @@ class StageContract:
269
  fail_blocks: bool = True
270
  run: Optional[Callable[[dict, dict], tuple[dict, dict]]] = None
271
  expectation: Optional[str] = None
 
272
 
273
  def resolve_rules(self, sample: dict) -> list[ThresholdRule]:
274
  if callable(self.rules):
@@ -302,6 +306,7 @@ def run_contract(contract: StageContract, sample: dict, stage_state: dict) -> St
302
  qc=qc,
303
  decision=qc.decision,
304
  data=data,
 
305
  )
306
 
307
 
 
86
  evaluate: Callable[[Any], QcStatus]
87
  expectation: Optional[str] = None
88
  optional: bool = False # missing metric -> WARN (not FAIL); used for conditional sub-checks
89
+ missing_detail: Optional[str] = None
90
 
91
  def apply(self, value: Any, detail: Optional[str] = None) -> Metric:
92
  return Metric(
 
154
  elif rule.optional:
155
  metrics.append(
156
  Metric(name=rule.metric, value=None, status=QcStatus.WARN, expected=rule.expectation,
157
+ detail=rule.missing_detail or "optional metric not evaluated (required evidence was not supplied)")
158
  )
159
  else:
160
  metrics.append(
 
228
  qc: Optional[QcResult] = None
229
  decision: Decision = Decision.CONTINUE
230
  data: dict = field(default_factory=dict) # stage-specific payload (parsed results)
231
+ evidence_level: str = "INTERNAL_COMPUTATION"
232
 
233
  def to_dict(self) -> dict:
234
  return {
 
243
  "qc": self.qc.to_dict() if self.qc else None,
244
  "decision": self.decision.value,
245
  "data": self.data,
246
+ "evidence_level": self.evidence_level,
247
  }
248
 
249
 
 
272
  fail_blocks: bool = True
273
  run: Optional[Callable[[dict, dict], tuple[dict, dict]]] = None
274
  expectation: Optional[str] = None
275
+ evidence_level: str = "INTERNAL_COMPUTATION"
276
 
277
  def resolve_rules(self, sample: dict) -> list[ThresholdRule]:
278
  if callable(self.rules):
 
306
  qc=qc,
307
  decision=qc.decision,
308
  data=data,
309
+ evidence_level=contract.evidence_level,
310
  )
311
 
312
 
app/ngs/execution.py ADDED
@@ -0,0 +1,265 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ """Real nf-core/sarek executor adapters.
2
+
3
+ No adapter falls back to the exploratory Python pipeline. Disabled or incomplete
4
+ infrastructure is reported as unavailable before a job is accepted.
5
+ """
6
+
7
+ from __future__ import annotations
8
+
9
+ import json
10
+ import os
11
+ import shutil
12
+ import subprocess
13
+ import sys
14
+ import uuid
15
+ from urllib.parse import urlparse
16
+ from datetime import datetime, timezone
17
+ from pathlib import Path
18
+ from typing import Any
19
+
20
+ from app.config import settings
21
+
22
+
23
+ def _now() -> str:
24
+ return datetime.now(timezone.utc).isoformat()
25
+
26
+
27
+ def _binary(name: str) -> bool:
28
+ return shutil.which(name) is not None
29
+
30
+
31
+ def executor_capabilities() -> dict[str, Any]:
32
+ local_missing = [name for name in ("nextflow", "java") if not _binary(name)]
33
+ if not (_binary("docker") or _binary("singularity") or _binary("apptainer")):
34
+ local_missing.append("docker, singularity, or apptainer")
35
+ slurm_missing = [name for name in ("nextflow", "java", "sbatch", "sacct") if not _binary(name)]
36
+ aws_missing = [
37
+ name for name, value in (
38
+ ("NGS_AWS_REGION", settings.NGS_AWS_REGION),
39
+ ("NGS_AWS_BATCH_JOB_QUEUE", settings.NGS_AWS_BATCH_JOB_QUEUE),
40
+ ("NGS_AWS_BATCH_JOB_DEFINITION", settings.NGS_AWS_BATCH_JOB_DEFINITION),
41
+ ) if not value
42
+ ]
43
+ try:
44
+ import boto3 # noqa: F401
45
+ except ImportError:
46
+ aws_missing.append("boto3")
47
+ return {
48
+ "workflow": {"name": "nf-core/sarek", "revision": "3.10.0"},
49
+ "executors": {
50
+ "local": {
51
+ "available": settings.NGS_LOCAL_EXECUTION_ENABLED and not local_missing,
52
+ "enabled": settings.NGS_LOCAL_EXECUTION_ENABLED,
53
+ "missing": local_missing,
54
+ "requirements": ["Nextflow", "Java", "Docker/Apptainer/Singularity", "durable local result directory"],
55
+ },
56
+ "slurm": {
57
+ "available": settings.NGS_SLURM_EXECUTION_ENABLED and not slurm_missing,
58
+ "enabled": settings.NGS_SLURM_EXECUTION_ENABLED,
59
+ "missing": slurm_missing,
60
+ "requirements": ["Nextflow head node", "sbatch/sacct", "shared samplesheet/reference/work/result paths"],
61
+ },
62
+ "awsbatch": {
63
+ "available": settings.NGS_AWS_BATCH_EXECUTION_ENABLED and not aws_missing,
64
+ "enabled": settings.NGS_AWS_BATCH_EXECUTION_ENABLED,
65
+ "missing": aws_missing,
66
+ "requirements": ["AWS Batch queue", "Nextflow driver job definition", "IAM role", "S3 inputs/results"],
67
+ },
68
+ },
69
+ "fallback": None,
70
+ "note": "Production executors never fall back to the exploratory preview.",
71
+ }
72
+
73
+
74
+ def _require(executor: str) -> None:
75
+ capability = executor_capabilities()["executors"][executor]
76
+ if not capability["available"]:
77
+ reasons = list(capability["missing"])
78
+ if not capability["enabled"]:
79
+ reasons.insert(0, f"NGS_{'AWS_BATCH' if executor == 'awsbatch' else executor.upper()}_EXECUTION_ENABLED")
80
+ raise RuntimeError("executor unavailable: " + ", ".join(reasons))
81
+
82
+
83
+ def _run_dir(run_id: str) -> Path:
84
+ path = Path(settings.NGS_RUN_ROOT).resolve() / run_id
85
+ path.mkdir(parents=True, exist_ok=False)
86
+ return path
87
+
88
+
89
+ def submit_run(executor: str, command_argv: list[str], outdir: str, user_id: str) -> dict[str, Any]:
90
+ if executor not in {"local", "slurm", "awsbatch"}:
91
+ raise ValueError(f"unsupported executor: {executor}")
92
+ _require(executor)
93
+ if executor == "local" and "-profile" in command_argv:
94
+ profile_index = command_argv.index("-profile") + 1
95
+ profile = command_argv[profile_index] if profile_index < len(command_argv) else ""
96
+ runtime = {"docker": "docker", "singularity": "singularity", "apptainer": "apptainer"}.get(profile)
97
+ if runtime and not _binary(runtime):
98
+ raise RuntimeError(f"executor unavailable: selected {runtime} runtime is not installed")
99
+ run_id = str(uuid.uuid4())
100
+ submitted_at = _now()
101
+
102
+ if executor == "local":
103
+ run_dir = _run_dir(run_id)
104
+ status_path = run_dir / "status.json"
105
+ log_path = run_dir / "nextflow.log"
106
+ worker_path = Path(__file__).with_name("nextflow_worker.py")
107
+ process = subprocess.Popen(
108
+ [sys.executable, str(worker_path), str(status_path), str(log_path), *command_argv],
109
+ cwd=str(run_dir), stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL, start_new_session=True,
110
+ )
111
+ executor_job_id = str(process.pid)
112
+ elif executor == "slurm":
113
+ run_dir = _run_dir(run_id)
114
+ script_path = run_dir / "launch.sh"
115
+ import shlex
116
+ script_path.write_text("#!/bin/sh\nset -eu\n" + shlex.join(command_argv) + "\n", encoding="utf-8")
117
+ script_path.chmod(0o700)
118
+ result = subprocess.run(["sbatch", "--parsable", str(script_path)], capture_output=True, text=True, timeout=30, check=True)
119
+ executor_job_id = result.stdout.strip().split(";")[0]
120
+ else:
121
+ import boto3
122
+ client = boto3.client("batch", region_name=settings.NGS_AWS_REGION)
123
+ response = client.submit_job(
124
+ jobName=f"bionexus-sarek-{run_id[:8]}",
125
+ jobQueue=settings.NGS_AWS_BATCH_JOB_QUEUE,
126
+ jobDefinition=settings.NGS_AWS_BATCH_JOB_DEFINITION,
127
+ containerOverrides={"command": command_argv, "environment": [{"name": "BIONEXUS_RUN_ID", "value": run_id}]},
128
+ tags={"BioNexusRunId": run_id, "Workflow": "nf-core-sarek-3.10.0"},
129
+ )
130
+ executor_job_id = response["jobId"]
131
+
132
+ record = {
133
+ "run_id": run_id, "state": "SUBMITTED", "executor": executor,
134
+ "executor_job_id": executor_job_id, "workflow": "nf-core/sarek", "revision": "3.10.0",
135
+ "outdir": outdir, "submitted_at": submitted_at, "updated_at": submitted_at,
136
+ "exit_code": None, "message": None, "user_id": user_id,
137
+ }
138
+ _persist_record(record)
139
+ return record
140
+
141
+
142
+ def _record_path(run_id: str) -> Path:
143
+ root = Path(settings.NGS_RUN_ROOT).resolve()
144
+ root.mkdir(parents=True, exist_ok=True)
145
+ return root / f"{run_id}.json"
146
+
147
+
148
+ def _persist_record(record: dict[str, Any]) -> None:
149
+ path = _record_path(record["run_id"])
150
+ temporary = path.with_suffix(".tmp")
151
+ temporary.write_text(json.dumps(record, sort_keys=True), encoding="utf-8")
152
+ temporary.replace(path)
153
+ try:
154
+ from app.services.supabase import get_supabase
155
+ get_supabase().table("ngs_production_runs").upsert(record, on_conflict="run_id").execute()
156
+ except Exception:
157
+ # The executor-local record keeps self-hosted operation functional; hosted
158
+ # deployments should apply migration 011 for durable cross-restart status.
159
+ pass
160
+
161
+
162
+ def _load_record(run_id: str, user_id: str) -> dict[str, Any]:
163
+ if not run_id or any(char not in "0123456789abcdef-" for char in run_id.lower()):
164
+ raise FileNotFoundError(run_id)
165
+ try:
166
+ from app.services.supabase import get_supabase
167
+ rows = get_supabase().table("ngs_production_runs").select("*").eq("run_id", run_id).eq("user_id", user_id).execute().data
168
+ if rows:
169
+ return dict(rows[0])
170
+ except Exception:
171
+ pass
172
+ path = _record_path(run_id)
173
+ if not path.is_file():
174
+ raise FileNotFoundError(run_id)
175
+ record = json.loads(path.read_text(encoding="utf-8"))
176
+ if record.get("user_id") != user_id:
177
+ raise FileNotFoundError(run_id)
178
+ return record
179
+
180
+
181
+ def get_run(run_id: str, user_id: str) -> dict[str, Any]:
182
+ record = _load_record(run_id, user_id)
183
+ executor = record["executor"]
184
+ try:
185
+ if executor == "local":
186
+ status_path = Path(settings.NGS_RUN_ROOT).resolve() / run_id / "status.json"
187
+ if status_path.is_file():
188
+ status = json.loads(status_path.read_text(encoding="utf-8"))
189
+ record.update({key: status.get(key) for key in ("state", "exit_code", "message")})
190
+ elif executor == "slurm":
191
+ result = subprocess.run(
192
+ ["sacct", "-j", record["executor_job_id"], "--noheader", "--parsable2", "--format=State,ExitCode"],
193
+ capture_output=True, text=True, timeout=20, check=True,
194
+ )
195
+ line = next((line for line in result.stdout.splitlines() if line.strip()), "")
196
+ raw_state, _, exit_code = line.partition("|")
197
+ record["state"] = _slurm_state(raw_state)
198
+ record["exit_code"] = int(exit_code.split(":")[0]) if exit_code else None
199
+ else:
200
+ import boto3
201
+ client = boto3.client("batch", region_name=settings.NGS_AWS_REGION)
202
+ jobs = client.describe_jobs(jobs=[record["executor_job_id"]]).get("jobs", [])
203
+ if jobs:
204
+ job = jobs[0]
205
+ record["state"] = _aws_state(job.get("status", "UNKNOWN"))
206
+ record["exit_code"] = job.get("container", {}).get("exitCode")
207
+ record["message"] = job.get("statusReason")
208
+ except Exception as exc:
209
+ record["message"] = f"Status refresh failed: {type(exc).__name__}"
210
+ record["updated_at"] = _now()
211
+ _persist_record(record)
212
+ return record
213
+
214
+
215
+ def _slurm_state(state: str) -> str:
216
+ state = state.upper().split()[0] if state else "UNKNOWN"
217
+ if state in {"PENDING", "CONFIGURING"}: return "PENDING"
218
+ if state in {"RUNNING", "COMPLETING"}: return "RUNNING"
219
+ if state == "COMPLETED": return "SUCCEEDED"
220
+ if state in {"FAILED", "CANCELLED", "TIMEOUT", "OUT_OF_MEMORY", "NODE_FAIL"}: return "FAILED"
221
+ return "UNKNOWN"
222
+
223
+
224
+ def _aws_state(state: str) -> str:
225
+ return {"SUBMITTED": "SUBMITTED", "PENDING": "PENDING", "RUNNABLE": "PENDING", "STARTING": "PENDING", "RUNNING": "RUNNING", "SUCCEEDED": "SUCCEEDED", "FAILED": "FAILED"}.get(state, "UNKNOWN")
226
+
227
+
228
+ def artifact_manifest(run_id: str, user_id: str) -> dict[str, Any]:
229
+ """Inventory actual output objects; never synthesize missing artifact evidence."""
230
+ record = _load_record(run_id, user_id)
231
+ outdir = record["outdir"]
232
+ if outdir.startswith("s3://"):
233
+ parsed = urlparse(outdir)
234
+ import boto3
235
+ client = boto3.client("s3", region_name=settings.NGS_AWS_REGION or None)
236
+ paginator = client.get_paginator("list_objects_v2")
237
+ files = [
238
+ f"s3://{parsed.netloc}/{item['Key']}"
239
+ for page in paginator.paginate(Bucket=parsed.netloc, Prefix=parsed.path.lstrip("/"))
240
+ for item in page.get("Contents", [])
241
+ ]
242
+ else:
243
+ root = Path(outdir).resolve()
244
+ if not root.is_dir():
245
+ files = []
246
+ else:
247
+ files = [str(path) for path in root.rglob("*") if path.is_file()]
248
+
249
+ lower = {path: path.lower() for path in files}
250
+ groups = {
251
+ "execution": [path for path, value in lower.items() if "pipeline_info/execution_" in value],
252
+ "multiqc": [path for path, value in lower.items() if "multiqc" in value],
253
+ "alignment": [path for path, value in lower.items() if value.endswith((".bam", ".cram", ".bai", ".crai"))],
254
+ "small_variants": [path for path, value in lower.items() if value.endswith((".vcf.gz", ".vcf.gz.tbi", ".g.vcf.gz", ".g.vcf.gz.tbi"))],
255
+ "coverage": [path for path, value in lower.items() if "mosdepth" in value or "coverage" in value],
256
+ "identity_qc": [path for path, value in lower.items() if any(token in value for token in ("contamination", "fingerprint", "verifybamid", "sex"))],
257
+ "provenance": [path for path, value in lower.items() if value.endswith(("run_manifest.json", "checksums.sha256"))],
258
+ }
259
+ missing = [name for name, evidence in groups.items() if not evidence]
260
+ return {
261
+ "run_id": run_id, "workflow": record["workflow"], "revision": record["revision"],
262
+ "source": outdir, "observed_file_count": len(files), "groups": groups,
263
+ "required_groups_complete": not missing, "missing_groups": missing,
264
+ "claim": "Observed artifact inventory only; presence does not imply QC or clinical validity.",
265
+ }
app/ngs/nextflow_worker.py ADDED
@@ -0,0 +1,33 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ """Detached local Nextflow process wrapper that leaves a durable status file."""
2
+
3
+ from __future__ import annotations
4
+
5
+ import json
6
+ import subprocess
7
+ import sys
8
+ from datetime import datetime, timezone
9
+ from pathlib import Path
10
+
11
+
12
+ def _write(path: Path, state: str, exit_code: int | None = None, message: str | None = None) -> None:
13
+ path.write_text(json.dumps({"state": state, "exit_code": exit_code, "message": message, "updated_at": datetime.now(timezone.utc).isoformat()}), encoding="utf-8")
14
+
15
+
16
+ def main() -> int:
17
+ if len(sys.argv) < 4:
18
+ return 2
19
+ status_path, log_path, command = Path(sys.argv[1]), Path(sys.argv[2]), sys.argv[3:]
20
+ _write(status_path, "RUNNING")
21
+ try:
22
+ with log_path.open("w", encoding="utf-8") as log:
23
+ completed = subprocess.run(command, stdout=log, stderr=subprocess.STDOUT, text=True, check=False)
24
+ state = "SUCCEEDED" if completed.returncode == 0 else "FAILED"
25
+ _write(status_path, state, completed.returncode, None if completed.returncode == 0 else "Nextflow exited non-zero; inspect nextflow.log")
26
+ return completed.returncode
27
+ except Exception as exc:
28
+ _write(status_path, "FAILED", None, f"{type(exc).__name__}: {exc}")
29
+ return 1
30
+
31
+
32
+ if __name__ == "__main__":
33
+ raise SystemExit(main())
app/ngs/orchestrator.py CHANGED
@@ -35,6 +35,58 @@ from app.ngs.contracts import (
35
 
36
  logger = logging.getLogger(__name__)
37
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
38
 
39
  class Pipeline:
40
  """Runs an ordered list of StageContracts over a shared sample/state context."""
@@ -69,10 +121,19 @@ class Pipeline:
69
  self.results.append(stage_result)
70
 
71
  if stage_result.qc:
72
- for m in stage_result.qc.metrics:
73
- if m.status == QcStatus.WARN:
 
 
 
 
 
 
 
 
 
74
  self.warnings.append(
75
- f"[{contract.step}] {m.name}: {m.value} ({m.expected})"
76
  )
77
 
78
  if stage_result.decision == Decision.STOP:
@@ -117,12 +178,16 @@ class Pipeline:
117
  "inputs": files,
118
  "reference": reference,
119
  "tools": [
120
- {"stage": r.step, "name": r.tool, "version": r.version}
 
121
  for r in self.results
122
  ],
123
  }
124
 
125
  def report(self) -> dict:
 
 
 
126
  return {
127
  "pipeline": self.name,
128
  "pipeline_version": self.version,
@@ -132,6 +197,16 @@ class Pipeline:
132
  "warnings": self.warnings,
133
  "stages": [r.to_dict() for r in self.results],
134
  "provenance": self.provenance,
 
 
 
 
 
 
 
 
 
 
135
  }
136
 
137
  def _overall_status(self) -> QcStatus:
 
35
 
36
  logger = logging.getLogger(__name__)
37
 
38
+ BENCHMARK_REGISTRY = [
39
+ {
40
+ "id": "giab-hg002-v4.2.1",
41
+ "source": "NIST Genome in a Bottle",
42
+ "scope": "HG002 germline small variants within benchmark regions",
43
+ "comparison_method": "GA4GH hap.py or vcfeval, stratified by variant type and genome context",
44
+ "url": "https://www.nist.gov/programs-projects/genome-bottle",
45
+ "status": "NOT_EVALUATED",
46
+ "metrics": None,
47
+ "reason": "This run did not supply HG002 reads, the matching truth VCF and benchmark BED, or a GA4GH comparison report.",
48
+ },
49
+ {
50
+ "id": "precisionfda-truth-v2",
51
+ "source": "FDA precisionFDA Truth Challenge V2",
52
+ "scope": "HG002/HG003/HG004 variants, including difficult-to-map regions",
53
+ "comparison_method": "Challenge-compatible precision, recall and F1 by region and variant class",
54
+ "url": "https://precision.fda.gov/challenges/10/results",
55
+ "status": "NOT_EVALUATED",
56
+ "metrics": None,
57
+ "reason": "No challenge-compatible callset and stratified evaluation output were produced by this run.",
58
+ },
59
+ {
60
+ "id": "seqc-gse47774",
61
+ "source": "NCBI GEO / SEQC consortium",
62
+ "scope": "RNA-seq accuracy and reproducibility reference dataset GSE47774",
63
+ "comparison_method": "Expression accuracy, replicate correlation and differential-expression reproducibility",
64
+ "url": "https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE47774",
65
+ "status": "NOT_APPLICABLE",
66
+ "metrics": None,
67
+ "reason": "This WGS/WES workflow does not produce RNA-seq expression estimates.",
68
+ },
69
+ ]
70
+
71
+ PRODUCTION_REQUIREMENTS = [
72
+ {"id": "complete_input", "label": "Complete FASTQ ingestion", "required": True,
73
+ "detail": "All records processed; no preview cap or silent subsampling."},
74
+ {"id": "production_alignment", "label": "Production read alignment", "required": True,
75
+ "detail": "Executed BWA-MEM2/DRAGMAP or another validated aligner with command and version."},
76
+ {"id": "bam_artifacts", "label": "Indexed alignment artifacts", "required": True,
77
+ "detail": "Coordinate-sorted BAM/CRAM, index, flagstat, idxstats and alignment metrics."},
78
+ {"id": "recalibration", "label": "Reference-matched recalibration", "required": True,
79
+ "detail": "BQSR or a documented no-BQSR workflow using build-matched known sites."},
80
+ {"id": "production_calls", "label": "Production germline callset", "required": True,
81
+ "detail": "Caller-generated VCF/gVCF with genotype, DP, GQ, AD, filters and index."},
82
+ {"id": "sample_qc", "label": "Sample QC and identity", "required": True,
83
+ "detail": "Coverage, duplication, insert size, contamination, sex and concordance when applicable."},
84
+ {"id": "truth_benchmark", "label": "GIAB/GA4GH truth comparison", "required": True,
85
+ "detail": "hap.py/vcfeval precision, recall and F1 inside the matching benchmark BED, stratified by SNP/INDEL."},
86
+ {"id": "reproducibility", "label": "Reproducible execution record", "required": True,
87
+ "detail": "Reference/resource checksums, exact commands, pipeline revision and container digests."},
88
+ ]
89
+
90
 
91
  class Pipeline:
92
  """Runs an ordered list of StageContracts over a shared sample/state context."""
 
121
  self.results.append(stage_result)
122
 
123
  if stage_result.qc:
124
+ unevaluated = stage_result.data.get("unevaluated")
125
+ if unevaluated:
126
+ self.warnings.append(
127
+ f"[{contract.step}] Not evaluated: {unevaluated}."
128
+ )
129
+ elif contract.step != "final_gate":
130
+ for metric in stage_result.qc.metrics:
131
+ if metric.status != QcStatus.WARN or metric.value is None:
132
+ continue
133
+ expected = f"; expected {metric.expected}" if metric.expected else ""
134
+ detail = f" — {metric.detail}" if metric.detail else ""
135
  self.warnings.append(
136
+ f"[{contract.step}] {metric.name}: {metric.value}{expected}{detail}"
137
  )
138
 
139
  if stage_result.decision == Decision.STOP:
 
178
  "inputs": files,
179
  "reference": reference,
180
  "tools": [
181
+ {"stage": r.step, "implementation": r.tool, "version": r.version,
182
+ "evidence_level": r.evidence_level}
183
  for r in self.results
184
  ],
185
  }
186
 
187
  def report(self) -> dict:
188
+ surrogate_stages = [r.step for r in self.results if r.evidence_level == "SURROGATE"]
189
+ requirements = [{**item, "status": "MISSING", "evidence": None}
190
+ for item in PRODUCTION_REQUIREMENTS]
191
  return {
192
  "pipeline": self.name,
193
  "pipeline_version": self.version,
 
197
  "warnings": self.warnings,
198
  "stages": [r.to_dict() for r in self.results],
199
  "provenance": self.provenance,
200
+ "validation": {
201
+ "claim": "NO_ACCURACY_CLAIM",
202
+ "analysis_grade": "EXPLORATORY_PREVIEW",
203
+ "research_ready": False,
204
+ "summary": "This internal sampled/surrogate workflow has not been validated against a public truth set.",
205
+ "same_or_better_supported": False,
206
+ "surrogate_stages": surrogate_stages,
207
+ "production_requirements": requirements,
208
+ "comparisons": BENCHMARK_REGISTRY,
209
+ },
210
  }
211
 
212
  def _overall_status(self) -> QcStatus:
app/ngs/portable_benchmark_report.json ADDED
@@ -0,0 +1,68 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ {
2
+ "benchmark": "bionexus-tiny-target-v1",
3
+ "classification": "SYNTHETIC_POSITIVE_CONTROL",
4
+ "status": "PASS",
5
+ "workflow_output_parity": true,
6
+ "biological_scope": "one heterozygous SNP in a 200 bp synthetic target",
7
+ "expected_call": {
8
+ "chrom": "chrTiny",
9
+ "pos": 50,
10
+ "ref": "C",
11
+ "alt": "G",
12
+ "genotype": "0/1",
13
+ "depth": 20,
14
+ "allelic_depth": "10,10"
15
+ },
16
+ "toolchain": {
17
+ "nextflow": "26.04.6",
18
+ "samtools": "1.24",
19
+ "bcftools": "1.24",
20
+ "filter": "QUAL>=20 && FORMAT/DP>=10"
21
+ },
22
+ "reports": [
23
+ {
24
+ "orchestrator": "Bio-Nexus direct",
25
+ "execution": "EXECUTED",
26
+ "normalized_sha256": "cefef322e336202da82549c1d5c09cf546a70b4aacfe9d1cc7e090a8aa23bbb1",
27
+ "tp": 1,
28
+ "fp": 0,
29
+ "fn": 0,
30
+ "precision": 1.0,
31
+ "recall": 1.0,
32
+ "f1": 1.0
33
+ },
34
+ {
35
+ "orchestrator": "Nextflow",
36
+ "execution": "EXECUTED",
37
+ "normalized_sha256": "cefef322e336202da82549c1d5c09cf546a70b4aacfe9d1cc7e090a8aa23bbb1",
38
+ "tp": 1,
39
+ "fp": 0,
40
+ "fn": 0,
41
+ "precision": 1.0,
42
+ "recall": 1.0,
43
+ "f1": 1.0
44
+ },
45
+ {
46
+ "orchestrator": "Galaxy wrapper command",
47
+ "execution": "EXECUTED_WITHOUT_GALAXY_SERVER",
48
+ "normalized_sha256": "cefef322e336202da82549c1d5c09cf546a70b4aacfe9d1cc7e090a8aa23bbb1",
49
+ "tp": 1,
50
+ "fp": 0,
51
+ "fn": 0,
52
+ "precision": 1.0,
53
+ "recall": 1.0,
54
+ "f1": 1.0
55
+ }
56
+ ],
57
+ "galaxy_validation": {
58
+ "wrapper_lint": "PASS",
59
+ "local_server_test": "NOT_COMPLETED",
60
+ "reason": "Galaxy 25.1 cloned, but its pinned dependency bootstrap was interrupted by proxy failures."
61
+ },
62
+ "limitations": [
63
+ "This validates deterministic workflow portability, not real-sample accuracy.",
64
+ "Only one synthetic heterozygous SNP in an easy 200 bp target is tested.",
65
+ "Indels, difficult regions, SVs, CNVs, contamination and annotation are not tested.",
66
+ "This does not establish parity with nf-core/sarek as a whole or a public Galaxy server."
67
+ ]
68
+ }
app/ngs/production.py ADDED
@@ -0,0 +1,189 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ """Production Sarek planning and fail-closed clinical software gates.
2
+
3
+ This module does not execute Nextflow inside the request process. It emits an auditable,
4
+ pinned launch contract for a durable compute worker and evaluates imported run evidence.
5
+ """
6
+
7
+ from __future__ import annotations
8
+
9
+ import hashlib
10
+ import hmac
11
+ import json
12
+ import re
13
+ import shlex
14
+ from typing import Any
15
+
16
+ from app.models.responses import NgsClinicalEvidenceRequest, NgsProductionPlanRequest
17
+ from app.config import settings
18
+
19
+ SAREK_RELEASE = "3.10.0"
20
+ SCHEMA_VERSION = "1.0"
21
+ SHA256_RE = re.compile(r"^[0-9a-fA-F]{64}$")
22
+
23
+
24
+ def _artifact(artifact_id: str, description: str, patterns: list[str]) -> dict[str, Any]:
25
+ return {"id": artifact_id, "required": True, "description": description, "patterns": patterns}
26
+
27
+
28
+ REQUIRED_ARTIFACTS = [
29
+ _artifact("execution", "Nextflow execution reports from the actual run", ["pipeline_info/execution_trace*.txt", "pipeline_info/execution_report*.html", "pipeline_info/execution_timeline*.html"]),
30
+ _artifact("multiqc", "Aggregated Sarek QC report and parsed data", ["multiqc/multiqc_report.html", "multiqc/multiqc_data/"]),
31
+ _artifact("alignment", "Coordinate-sorted alignment and index", ["**/*.bam|**/*.cram", "**/*.bai|**/*.crai"]),
32
+ _artifact("small_variants", "Caller VCF/gVCF and tabix index", ["**/*.vcf.gz", "**/*.vcf.gz.tbi"]),
33
+ _artifact("coverage", "Assay-appropriate coverage metrics", ["reports/**/mosdepth*", "reports/**/coverage*"]),
34
+ _artifact("identity_qc", "Sample identity, contamination and sex/ploidy evidence", ["reports/**/*contamination*", "reports/**/*fingerprint*", "reports/**/*sex*"]),
35
+ _artifact("provenance", "Checksums and immutable workflow/resource identifiers", ["run_manifest.json", "checksums.sha256"]),
36
+ ]
37
+
38
+ PROVENANCE_REQUIREMENTS = [
39
+ "nf-core/sarek revision and Nextflow version",
40
+ "container image digest for every executed process",
41
+ "input sample-sheet SHA-256 and input file checksums",
42
+ "reference FASTA, indexes, known-sites, target BED and annotation-cache checksums",
43
+ "complete command argument array, profile and custom configuration checksum",
44
+ "execution trace with process exit status and retry history",
45
+ ]
46
+
47
+
48
+ def build_production_plan(request: NgsProductionPlanRequest) -> dict[str, Any]:
49
+ blockers: list[str] = []
50
+ warnings: list[str] = []
51
+ if request.assay == "WES" and not (request.target_bed or "").strip():
52
+ blockers.append("WES requires a reference-matched target BED.")
53
+ if request.execution_profile in {"slurm", "awsbatch"} and not (request.custom_config or "").strip():
54
+ blockers.append(f"{request.execution_profile} requires a reviewed Nextflow custom configuration.")
55
+ if request.input_type == "FASTQ" and request.start_step != "mapping":
56
+ blockers.append("FASTQ input must start at Sarek's mapping step.")
57
+ if request.input_type == "CRAM" and request.start_step == "mapping":
58
+ blockers.append("CRAM input must start at markduplicates or variant_calling and include its CRAI in the sample sheet.")
59
+ if request.input_type == "BAM" and request.start_step == "mapping":
60
+ warnings.append("BAM at the mapping step must be an unmapped BAM; prepared alignments should start at markduplicates or variant_calling with an index.")
61
+ if request.genome == "GRCh37":
62
+ warnings.append("GRCh37 resources must all use the same contig naming and bundle; hg19 resources are not interchangeable by label alone.")
63
+ if request.clinical_intent:
64
+ warnings.append("Clinical intent activates a fail-closed evidence gate; this launch plan alone cannot validate or authorize an assay.")
65
+ if request.sample_model != "singleton" and request.caller != "haplotypecaller":
66
+ warnings.append("The selected caller produces per-sample results here; joint germline calling is enabled only for HaplotypeCaller.")
67
+
68
+ profile = request.execution_profile
69
+ sarek_genome = {"GRCh38": "GATK.GRCh38", "GRCh37": "GATK.GRCh37"}[request.genome]
70
+ argv = [
71
+ "nextflow", "run", "nf-core/sarek", "-r", SAREK_RELEASE,
72
+ "-profile", profile,
73
+ "--input", request.samplesheet_path,
74
+ "--outdir", request.outdir,
75
+ "--genome", sarek_genome,
76
+ "--tools", f"{request.caller}{',vep' if request.annotate_with_vep else ''}",
77
+ "--step", request.start_step,
78
+ ]
79
+ if request.assay == "WES":
80
+ argv.extend(["--wes", "--intervals", request.target_bed or ""])
81
+ if request.custom_config:
82
+ argv.extend(["-c", request.custom_config])
83
+ if request.sample_model != "singleton" and request.caller == "haplotypecaller":
84
+ argv.append("--joint_germline")
85
+
86
+ ready = not blockers
87
+ return {
88
+ "schema_version": SCHEMA_VERSION,
89
+ "workflow": {
90
+ "name": "nf-core/sarek",
91
+ "revision": SAREK_RELEASE,
92
+ "execution_engine": "Nextflow",
93
+ "assay": request.assay,
94
+ "sample_model": request.sample_model,
95
+ "input_type": request.input_type,
96
+ "start_step": request.start_step,
97
+ },
98
+ "state": "PLANNED" if ready else "BLOCKED",
99
+ "ready_to_launch": ready,
100
+ "blockers": blockers,
101
+ "warnings": warnings,
102
+ "command_argv": argv,
103
+ "command_display": shlex.join(argv),
104
+ "required_artifacts": REQUIRED_ARTIFACTS,
105
+ "provenance_requirements": PROVENANCE_REQUIREMENTS,
106
+ "clinical_boundary": {
107
+ "clinical_intent": request.clinical_intent,
108
+ "current_status": "NOT_CLINICALLY_RELEASABLE",
109
+ "reason": "No executed run evidence, assay-validation record, benchmark evidence or authorized review has been imported.",
110
+ },
111
+ }
112
+
113
+
114
+ def _present(value: str | None) -> bool:
115
+ return bool(value and value.strip())
116
+
117
+
118
+ def _sha256(value: str | None) -> bool:
119
+ return bool(value and SHA256_RE.fullmatch(value))
120
+
121
+
122
+ def clinical_signature_payload(evidence: NgsClinicalEvidenceRequest) -> bytes:
123
+ """Canonical bytes signed by the trusted evidence-import service."""
124
+ if hasattr(evidence, "model_dump"):
125
+ payload = evidence.model_dump(mode="json", exclude={"evidence_signature"})
126
+ else: # Pydantic v1 compatibility for older self-hosted deployments.
127
+ payload = evidence.dict(exclude={"evidence_signature"})
128
+ return json.dumps(payload, sort_keys=True, separators=(",", ":")).encode("utf-8")
129
+
130
+
131
+ def _signature_verified(evidence: NgsClinicalEvidenceRequest) -> bool:
132
+ key = settings.NGS_CLINICAL_EVIDENCE_HMAC_KEY
133
+ if not key or not evidence.evidence_signature or not _sha256(evidence.evidence_bundle_sha256):
134
+ return False
135
+ expected = hmac.new(key.encode("utf-8"), clinical_signature_payload(evidence), hashlib.sha256).hexdigest()
136
+ return hmac.compare_digest(expected, evidence.evidence_signature)
137
+
138
+
139
+ def evaluate_clinical_evidence(evidence: NgsClinicalEvidenceRequest) -> dict[str, Any]:
140
+ metrics_present = all(value is not None for value in (
141
+ evidence.snv_precision, evidence.snv_recall,
142
+ evidence.indel_precision, evidence.indel_recall,
143
+ ))
144
+ truth_is_external = _present(evidence.truthset_name) and "synthetic" not in evidence.truthset_name.lower()
145
+ review_complete = evidence.human_reviewed and _present(evidence.reviewer_id) and _present(evidence.release_signature_id)
146
+ signature_verified = _signature_verified(evidence)
147
+
148
+ checks: list[tuple[str, str, bool, str]] = [
149
+ ("evidence_signature", "Trusted evidence bundle signature verified", signature_verified, evidence.evidence_bundle_sha256 or "missing or unverified"),
150
+ ("assay_validation", "Validated assay record", _present(evidence.assay_validation_id), evidence.assay_validation_id or "missing"),
151
+ ("execution", "Pinned workflow completed", evidence.workflow_status == "COMPLETED" and _present(evidence.sarek_revision), f"{evidence.workflow_status}; revision={evidence.sarek_revision or 'missing'}"),
152
+ ("reference", "Reference and sample-sheet checksums", evidence.reference_build is not None and _sha256(evidence.reference_manifest_sha256) and _sha256(evidence.samplesheet_sha256), evidence.reference_build or "missing"),
153
+ ("containers", "Container digests captured", evidence.container_digests_complete, str(evidence.container_digests_complete)),
154
+ ("complete_input", "All input records processed", evidence.complete_input_processed, str(evidence.complete_input_processed)),
155
+ ("artifacts", "Required artifacts imported", evidence.required_artifacts_present, str(evidence.required_artifacts_present)),
156
+ ("qc", "Run QC passed", evidence.qc_pass, str(evidence.qc_pass)),
157
+ ("identity", "Sample identity passed", evidence.sample_identity_pass, str(evidence.sample_identity_pass)),
158
+ ("contamination", "Contamination gate passed", evidence.contamination_pass, str(evidence.contamination_pass)),
159
+ ("sex_ploidy", "Sex/ploidy reviewed", evidence.sex_ploidy_reviewed, str(evidence.sex_ploidy_reviewed)),
160
+ ("truth_benchmark", "External truth benchmark passed its approved protocol", bool(truth_is_external and _present(evidence.benchmark_protocol_id) and evidence.benchmark_acceptance_pass and _sha256(evidence.confident_regions_sha256) and evidence.same_sample_reference_regions and metrics_present), evidence.truthset_name or "missing"),
161
+ ("human_release", "Authorized human review and release signature", review_complete, evidence.reviewer_id or "missing"),
162
+ ("deviations", "No unresolved deviations", not evidence.unresolved_deviations, ", ".join(evidence.unresolved_deviations) or "none"),
163
+ ]
164
+ gates = [
165
+ {"id": gate_id, "label": label, "status": "PASS" if passed else "FAIL", "evidence": detail}
166
+ for gate_id, label, passed, detail in checks
167
+ ]
168
+ missing = [label for _gate_id, label, passed, _detail in checks if not passed]
169
+ passed = not missing
170
+ return {
171
+ "schema_version": SCHEMA_VERSION,
172
+ "status": "SOFTWARE_GATE_PASSED" if passed else "NOT_CLINICALLY_RELEASABLE",
173
+ "clinically_validated": False,
174
+ "gates": gates,
175
+ "missing_or_failed": missing,
176
+ "benchmark_summary": {
177
+ "truthset": evidence.truthset_name,
178
+ "protocol_id": evidence.benchmark_protocol_id,
179
+ "same_sample_reference_regions": evidence.same_sample_reference_regions,
180
+ "acceptance_pass": evidence.benchmark_acceptance_pass,
181
+ "metrics": {
182
+ "snv_precision": evidence.snv_precision,
183
+ "snv_recall": evidence.snv_recall,
184
+ "indel_precision": evidence.indel_precision,
185
+ "indel_recall": evidence.indel_recall,
186
+ },
187
+ },
188
+ "disclaimer": "This is a software evidence gate, not clinical validation, accreditation, diagnosis or report authorization. Laboratory validation and jurisdiction-specific governance remain required.",
189
+ }
app/ngs/stages/stage10_identity.py CHANGED
@@ -107,7 +107,7 @@ def _stage10_run(sample: dict, state: dict) -> tuple[dict, dict]:
107
  if records is None:
108
  return ({"error": "identity needs aligned reads",
109
  "unevaluated": "no aligned reads"},
110
- {"identity_checked": 0.0}) # WARN: cannot assert identity without alignment
111
  eng = identity_engine(
112
  records,
113
  expected_gt=sample.get("expected_gt"),
@@ -120,16 +120,37 @@ def _stage10_run(sample: dict, state: dict) -> tuple[dict, dict]:
120
  # unevaluated (WARN) rather than silently passing or hard-failing the pipeline.
121
  if eng["concordance"] is None and eng["sex_match"] is None:
122
  eng["unevaluated"] = ("no expected genotype or sex supplied; sample identity not asserted")
123
- return eng, {"identity_checked": 0.0}
 
124
  metrics = {}
125
  if eng["concordance"] is not None:
126
  metrics["concordance_ok"] = eng["concordance"]
127
  if eng["sex_match"] is not None:
128
  metrics["sex_ok"] = 100.0 if eng["sex_match"] else 0.0
129
- metrics["identity_checked"] = 100.0
130
  return eng, metrics
131
 
132
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
133
  def stage10_contract() -> StageContract:
134
  return StageContract(
135
  step="identity",
@@ -137,24 +158,12 @@ def stage10_contract() -> StageContract:
137
  version="0.1.0",
138
  inputs=["processed_bam", "expected_genotypes?", "expected_sex?"],
139
  outputs=["identity_report"],
140
- rules=[
141
- ThresholdRule(name="identity_checked", metric="identity_checked",
142
- evaluate=lambda v: _checked_rule(v)),
143
- ThresholdRule(name="concordance_ok", metric="concordance_ok",
144
- evaluate=lambda v: _match_rule(v, 99.0), optional=True),
145
- ThresholdRule(name="sex_ok", metric="sex_ok",
146
- evaluate=lambda v: _bool_rule(v), optional=True),
147
- ],
148
  fail_blocks=True, # mislabeled/swapped sample stops the run
149
  run=_stage10_run,
150
  )
151
 
152
 
153
- def _checked_rule(v):
154
- # 100 = a real check ran; 0 = unevaluated (WARN, not blocking)
155
- return QcStatus.PASS if float(v) >= 100.0 else QcStatus.WARN
156
-
157
-
158
  def _match_rule(v, ok):
159
  try:
160
  v = float(v)
@@ -177,16 +186,13 @@ def run_identity(
177
  from app.ngs.contracts import apply_rules, QcResult
178
  eng = identity_engine(records, expected_gt, chr_x_depth, chr_y_depth, expected_sex)
179
  metrics = {}
180
- if eng["concordance"] is None and eng["sex_match"] is None:
181
- metrics["identity_checked"] = 0.0
182
- else:
183
- metrics["identity_checked"] = 100.0
184
  if eng["concordance"] is not None:
185
  metrics["concordance_ok"] = eng["concordance"]
186
  if eng["sex_match"] is not None:
187
  metrics["sex_ok"] = 100.0 if eng["sex_match"] else 0.0
188
  contract = stage10_contract()
189
- result = QcResult.from_metrics(apply_rules(contract.resolve_rules({}), metrics),
 
190
  fail_blocks=True)
191
  return {
192
  "result": {"step": "identity", "qc": result.to_dict(),
 
107
  if records is None:
108
  return ({"error": "identity needs aligned reads",
109
  "unevaluated": "no aligned reads"},
110
+ {})
111
  eng = identity_engine(
112
  records,
113
  expected_gt=sample.get("expected_gt"),
 
120
  # unevaluated (WARN) rather than silently passing or hard-failing the pipeline.
121
  if eng["concordance"] is None and eng["sex_match"] is None:
122
  eng["unevaluated"] = ("no expected genotype or sex supplied; sample identity not asserted")
123
+ eng["status"] = "NOT_EVALUATED"
124
+ return eng, {}
125
  metrics = {}
126
  if eng["concordance"] is not None:
127
  metrics["concordance_ok"] = eng["concordance"]
128
  if eng["sex_match"] is not None:
129
  metrics["sex_ok"] = 100.0 if eng["sex_match"] else 0.0
 
130
  return eng, metrics
131
 
132
 
133
+ def _identity_rules(sample: dict) -> list[ThresholdRule]:
134
+ rules: list[ThresholdRule] = []
135
+ if sample.get("expected_gt"):
136
+ rules.append(ThresholdRule(
137
+ name="concordance_ok", metric="concordance_ok",
138
+ evaluate=lambda v: _match_rule(v, 99.0), expectation=">= 99%",
139
+ ))
140
+ if sample.get("expected_sex"):
141
+ rules.append(ThresholdRule(
142
+ name="sex_ok", metric="sex_ok", evaluate=lambda v: _bool_rule(v),
143
+ expectation="matches declared sample sex",
144
+ ))
145
+ if not rules:
146
+ rules.append(ThresholdRule(
147
+ name="identity_evidence", metric="identity_evidence",
148
+ evaluate=lambda _v: QcStatus.PASS, optional=True,
149
+ missing_detail="Not evaluated: expected genotypes and declared sample sex were not supplied.",
150
+ ))
151
+ return rules
152
+
153
+
154
  def stage10_contract() -> StageContract:
155
  return StageContract(
156
  step="identity",
 
158
  version="0.1.0",
159
  inputs=["processed_bam", "expected_genotypes?", "expected_sex?"],
160
  outputs=["identity_report"],
161
+ rules=_identity_rules,
 
 
 
 
 
 
 
162
  fail_blocks=True, # mislabeled/swapped sample stops the run
163
  run=_stage10_run,
164
  )
165
 
166
 
 
 
 
 
 
167
  def _match_rule(v, ok):
168
  try:
169
  v = float(v)
 
186
  from app.ngs.contracts import apply_rules, QcResult
187
  eng = identity_engine(records, expected_gt, chr_x_depth, chr_y_depth, expected_sex)
188
  metrics = {}
 
 
 
 
189
  if eng["concordance"] is not None:
190
  metrics["concordance_ok"] = eng["concordance"]
191
  if eng["sex_match"] is not None:
192
  metrics["sex_ok"] = 100.0 if eng["sex_match"] else 0.0
193
  contract = stage10_contract()
194
+ sample = {"expected_gt": expected_gt, "expected_sex": expected_sex}
195
+ result = QcResult.from_metrics(apply_rules(contract.resolve_rules(sample), metrics),
196
  fail_blocks=True)
197
  return {
198
  "result": {"step": "identity", "qc": result.to_dict(),
app/ngs/stages/stage11_variant_calling.py CHANGED
@@ -131,6 +131,7 @@ def stage11_contract() -> StageContract:
131
  rules=[],
132
  fail_blocks=False,
133
  run=_stage11_run,
 
134
  )
135
 
136
 
 
131
  rules=[],
132
  fail_blocks=False,
133
  run=_stage11_run,
134
+ evidence_level="SURROGATE",
135
  )
136
 
137
 
app/ngs/stages/stage14_filter.py CHANGED
@@ -63,15 +63,14 @@ def filter_variants(variants: list[dict], max_af: float = 0.01) -> dict:
63
 
64
 
65
  def _stage14_run(sample: dict, state: dict) -> tuple[dict, dict]:
66
- variants = state.get("variants", {}).get("qc", {}).get("variants") or \
67
- state.get("variants", {}).get("normalized")
68
  if variants is None:
69
- return {"error": "filtering needs QC'd variants"}, {"candidate_frac": 0.0}
 
 
70
  filtered = filter_variants(variants)
71
  state.setdefault("variants", {})["filtered"] = filtered
72
- total = filtered["n_final"] + filtered["n_rejected"]
73
- frac = filtered["n_final"] / total * 100.0 if total else 0.0
74
- return filtered, {"candidate_frac": frac}
75
 
76
 
77
  def stage14_contract() -> StageContract:
@@ -82,31 +81,29 @@ def stage14_contract() -> StageContract:
82
  inputs=["qc_variants"],
83
  outputs=["final_variants"],
84
  rules=[
85
- ThresholdRule(name="candidate_frac", metric="candidate_frac",
86
- evaluate=lambda v: _range_rule(v, 0.1, 30.0)),
87
  ],
88
  fail_blocks=False,
89
  run=_stage14_run,
90
  )
91
 
92
 
93
- def _range_rule(v, lo, hi):
94
  from app.ngs.contracts import QcStatus
95
  try:
96
  v = float(v)
97
  except (TypeError, ValueError):
98
  return QcStatus.FAIL
99
- return QcStatus.PASS if lo <= v <= hi else QcStatus.WARN
100
 
101
 
102
  def run_variant_filter(variants: list[dict], max_af: float = 0.01) -> dict:
103
  from app.ngs.contracts import apply_rules, QcResult
104
  filtered = filter_variants(variants, max_af=max_af)
105
- total = filtered["n_final"] + filtered["n_rejected"]
106
- frac = filtered["n_final"] / total * 100.0 if total else 0.0
107
  contract = stage14_contract()
108
  result = QcResult.from_metrics(
109
- apply_rules(contract.resolve_rules({}), {"candidate_frac": frac}), fail_blocks=False)
110
  return {
111
  "result": {"step": "variant_filter", "qc": result.to_dict(),
112
  "decision": result.decision.value,
 
63
 
64
 
65
  def _stage14_run(sample: dict, state: dict) -> tuple[dict, dict]:
66
+ variants = state.get("variants", {}).get("qc", {}).get("variants")
 
67
  if variants is None:
68
+ variants = state.get("variants", {}).get("normalized")
69
+ if variants is None:
70
+ return {"error": "filtering needs QC'd variants"}, {"filter_completed": 0.0}
71
  filtered = filter_variants(variants)
72
  state.setdefault("variants", {})["filtered"] = filtered
73
+ return filtered, {"filter_completed": 100.0}
 
 
74
 
75
 
76
  def stage14_contract() -> StageContract:
 
81
  inputs=["qc_variants"],
82
  outputs=["final_variants"],
83
  rules=[
84
+ ThresholdRule(name="filter_completed", metric="filter_completed",
85
+ evaluate=lambda v: _completed_rule(v), expectation="completed"),
86
  ],
87
  fail_blocks=False,
88
  run=_stage14_run,
89
  )
90
 
91
 
92
+ def _completed_rule(v):
93
  from app.ngs.contracts import QcStatus
94
  try:
95
  v = float(v)
96
  except (TypeError, ValueError):
97
  return QcStatus.FAIL
98
+ return QcStatus.PASS if v >= 100.0 else QcStatus.FAIL
99
 
100
 
101
  def run_variant_filter(variants: list[dict], max_af: float = 0.01) -> dict:
102
  from app.ngs.contracts import apply_rules, QcResult
103
  filtered = filter_variants(variants, max_af=max_af)
 
 
104
  contract = stage14_contract()
105
  result = QcResult.from_metrics(
106
+ apply_rules(contract.resolve_rules({}), {"filter_completed": 100.0}), fail_blocks=False)
107
  return {
108
  "result": {"step": "variant_filter", "qc": result.to_dict(),
109
  "decision": result.decision.value,
app/ngs/stages/stage15_sv.py CHANGED
@@ -123,6 +123,7 @@ def stage15_contract() -> StageContract:
123
  ],
124
  fail_blocks=False,
125
  run=_stage15_run,
 
126
  )
127
 
128
 
 
123
  ],
124
  fail_blocks=False,
125
  run=_stage15_run,
126
+ evidence_level="SURROGATE",
127
  )
128
 
129
 
app/ngs/stages/stage16_cnv.py CHANGED
@@ -132,6 +132,7 @@ def stage16_contract() -> StageContract:
132
  ],
133
  fail_blocks=False,
134
  run=_stage16_run,
 
135
  )
136
 
137
 
 
132
  ],
133
  fail_blocks=False,
134
  run=_stage16_run,
135
+ evidence_level="SURROGATE",
136
  )
137
 
138
 
app/ngs/stages/stage17_annotation.py CHANGED
@@ -155,15 +155,17 @@ def run_annotation(variants: list[dict], transcripts: list[dict]) -> dict:
155
 
156
 
157
  def _stage17_run(sample: dict, state: dict) -> tuple[dict, dict]:
158
- variants = state.get("variants", {}).get("filtered", {}).get("final") or \
159
- state.get("variants", {}).get("final")
 
160
  transcripts = sample.get("transcripts") or state.get("transcripts") or []
161
  if variants is None:
162
- return {"error": "annotation needs final variants"}, {"annotate_frac": 0.0}
163
  report = run_annotation(variants, transcripts)
 
 
164
  state.setdefault("variants", {})["annotated"] = report
165
- frac = report["n_annotated"] / len(variants) * 100.0 if variants else 100.0
166
- return report, {"annotate_frac": frac}
167
 
168
 
169
  def stage17_contract() -> StageContract:
@@ -174,8 +176,8 @@ def stage17_contract() -> StageContract:
174
  inputs=["final_variants", "transcripts"],
175
  outputs=["annotated_variants"],
176
  rules=[
177
- ThresholdRule(name="annotate_frac", metric="annotate_frac",
178
- evaluate=lambda v: _pct_rule(v, 90, 70)),
179
  ],
180
  fail_blocks=False,
181
  run=_stage17_run,
@@ -197,10 +199,9 @@ def _pct_rule(v, ok, warn):
197
  def run_annotation_stage(variants: list[dict], transcripts: list[dict]) -> dict:
198
  from app.ngs.contracts import apply_rules, QcResult
199
  report = run_annotation(variants, transcripts)
200
- frac = report["n_annotated"] / len(variants) * 100.0 if variants else 100.0
201
  contract = stage17_contract()
202
  result = QcResult.from_metrics(
203
- apply_rules(contract.resolve_rules({}), {"annotate_frac": frac}), fail_blocks=False)
204
  return {
205
  "result": {"step": "annotation", "qc": result.to_dict(),
206
  "decision": result.decision.value, "data": report["summary"]},
 
155
 
156
 
157
  def _stage17_run(sample: dict, state: dict) -> tuple[dict, dict]:
158
+ variants = state.get("variants", {}).get("filtered", {}).get("final")
159
+ if variants is None:
160
+ variants = state.get("variants", {}).get("final")
161
  transcripts = sample.get("transcripts") or state.get("transcripts") or []
162
  if variants is None:
163
+ return {"error": "annotation needs final variants"}, {"annotation_completed": 0.0}
164
  report = run_annotation(variants, transcripts)
165
+ report["status"] = "NOT_APPLICABLE" if not variants else "COMPLETED"
166
+ report["n_input_variants"] = len(variants)
167
  state.setdefault("variants", {})["annotated"] = report
168
+ return report, {"annotation_completed": 100.0}
 
169
 
170
 
171
  def stage17_contract() -> StageContract:
 
176
  inputs=["final_variants", "transcripts"],
177
  outputs=["annotated_variants"],
178
  rules=[
179
+ ThresholdRule(name="annotation_completed", metric="annotation_completed",
180
+ evaluate=lambda v: _pct_rule(v, 100, 100), expectation="completed"),
181
  ],
182
  fail_blocks=False,
183
  run=_stage17_run,
 
199
  def run_annotation_stage(variants: list[dict], transcripts: list[dict]) -> dict:
200
  from app.ngs.contracts import apply_rules, QcResult
201
  report = run_annotation(variants, transcripts)
 
202
  contract = stage17_contract()
203
  result = QcResult.from_metrics(
204
+ apply_rules(contract.resolve_rules({}), {"annotation_completed": 100.0}), fail_blocks=False)
205
  return {
206
  "result": {"step": "annotation", "qc": result.to_dict(),
207
  "decision": result.decision.value, "data": report["summary"]},
app/ngs/stages/stage21_final_gate.py CHANGED
@@ -75,19 +75,17 @@ def evaluate_gate(pipeline_report: dict) -> dict:
75
  f"{verdict}: {len(blocked)} blocking gate(s) failed"
76
  if blocked else
77
  (f"{verdict}: proceeding with {len(warners)} warning(s)"
78
- if warners else "all {passed} pipeline stages passed")),
79
  }
80
 
81
 
82
  def _stage21_run(sample: dict, state: dict) -> tuple[dict, dict]:
83
  report = sample.get("pipeline_report") or state.get("pipeline_report")
84
  if report is None:
85
- return {"error": "final gate needs the pipeline report"}, {"gate_ok": 100.0}
86
  gate = evaluate_gate(report)
87
  state.setdefault("final_gate", {})["result"] = gate
88
- ok = 100.0 if gate["verdict"] == "ANALYSIS_READY" else (
89
- 60.0 if gate["verdict"] == "ANALYSIS_READY_WITH_WARNINGS" else 0.0)
90
- return gate, {"gate_ok": ok}
91
 
92
 
93
  def stage21_contract() -> StageContract:
@@ -98,8 +96,11 @@ def stage21_contract() -> StageContract:
98
  inputs=["pipeline_report"],
99
  outputs=["readiness_verdict"],
100
  rules=[
101
- ThresholdRule(name="gate_ok", metric="gate_ok",
102
- evaluate=lambda v: _gate_rule(v)),
 
 
 
103
  ],
104
  fail_blocks=True, # a NOT_READY verdict stops the report from being finalised
105
  run=_stage21_run,
@@ -107,13 +108,9 @@ def stage21_contract() -> StageContract:
107
 
108
 
109
  def _gate_rule(v):
110
- try:
111
- v = float(v)
112
- except (TypeError, ValueError):
113
- return QcStatus.FAIL
114
- if v >= 100.0:
115
  return QcStatus.PASS
116
- if v >= 60.0:
117
  return QcStatus.WARN
118
  return QcStatus.FAIL
119
 
@@ -121,10 +118,9 @@ def _gate_rule(v):
121
  def run_final_gate(pipeline_report: dict) -> dict:
122
  from app.ngs.contracts import apply_rules, QcResult
123
  gate = evaluate_gate(pipeline_report)
124
- ok = 100.0 if gate["verdict"] == "ANALYSIS_READY" else (
125
- 60.0 if gate["verdict"] == "ANALYSIS_READY_WITH_WARNINGS" else 0.0)
126
  contract = stage21_contract()
127
- result = QcResult.from_metrics(apply_rules(contract.resolve_rules({}), {"gate_ok": ok}),
 
128
  fail_blocks=True)
129
  return {
130
  "result": {"step": "final_gate", "qc": result.to_dict(),
 
75
  f"{verdict}: {len(blocked)} blocking gate(s) failed"
76
  if blocked else
77
  (f"{verdict}: proceeding with {len(warners)} warning(s)"
78
+ if warners else f"all {passed} pipeline stages passed")),
79
  }
80
 
81
 
82
  def _stage21_run(sample: dict, state: dict) -> tuple[dict, dict]:
83
  report = sample.get("pipeline_report") or state.get("pipeline_report")
84
  if report is None:
85
+ return {"error": "final gate needs the pipeline report"}, {"readiness_verdict": "NOT_ANALYSIS_READY"}
86
  gate = evaluate_gate(report)
87
  state.setdefault("final_gate", {})["result"] = gate
88
+ return gate, {"readiness_verdict": gate["verdict"]}
 
 
89
 
90
 
91
  def stage21_contract() -> StageContract:
 
96
  inputs=["pipeline_report"],
97
  outputs=["readiness_verdict"],
98
  rules=[
99
+ ThresholdRule(
100
+ name="readiness_verdict", metric="readiness_verdict",
101
+ evaluate=lambda v: _gate_rule(v),
102
+ expectation="ANALYSIS_READY or ANALYSIS_READY_WITH_WARNINGS",
103
+ ),
104
  ],
105
  fail_blocks=True, # a NOT_READY verdict stops the report from being finalised
106
  run=_stage21_run,
 
108
 
109
 
110
  def _gate_rule(v):
111
+ if v == "ANALYSIS_READY":
 
 
 
 
112
  return QcStatus.PASS
113
+ if v == "ANALYSIS_READY_WITH_WARNINGS":
114
  return QcStatus.WARN
115
  return QcStatus.FAIL
116
 
 
118
  def run_final_gate(pipeline_report: dict) -> dict:
119
  from app.ngs.contracts import apply_rules, QcResult
120
  gate = evaluate_gate(pipeline_report)
 
 
121
  contract = stage21_contract()
122
+ result = QcResult.from_metrics(apply_rules(
123
+ contract.resolve_rules({}), {"readiness_verdict": gate["verdict"]}),
124
  fail_blocks=True)
125
  return {
126
  "result": {"step": "final_gate", "qc": result.to_dict(),
app/ngs/stages/stage4_reference.py CHANGED
@@ -40,6 +40,23 @@ def _stage4_run(sample: dict, state: dict) -> tuple[dict, dict]:
40
  meta = sample.get("metadata") or {}
41
  annotation_build = sample.get("annotation_build") or meta.get("annotation_build")
42
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
43
  if not ref_id:
44
  return {"error": "no reference specified"}, {"reference_resolved": 0.0,
45
  "build_compatible": 100.0}
 
40
  meta = sample.get("metadata") or {}
41
  annotation_build = sample.get("annotation_build") or meta.get("annotation_build")
42
 
43
+ if sample.get("synthetic_reference"):
44
+ ref = {
45
+ "id": "synthetic-positive-control",
46
+ "reference_kind": "SYNTHETIC",
47
+ "build": "NOT_APPLICABLE",
48
+ "requested_template": str(ref_id or "not specified"),
49
+ "description": "Deterministic reference generated from demonstration reads",
50
+ "species": "NOT_APPLICABLE",
51
+ "generated_for_this_run": True,
52
+ "bundled": False,
53
+ "artifacts": [{"name": "In-memory synthetic reference", "present": True}],
54
+ "build_compatible": "NOT_APPLICABLE",
55
+ "build_message": "Synthetic coordinates are not GRCh37/GRCh38 coordinates.",
56
+ }
57
+ state.setdefault("reference", {})["declared"] = ref
58
+ return ref, {"reference_resolved": 100.0, "build_compatible": 100.0}
59
+
60
  if not ref_id:
61
  return {"error": "no reference specified"}, {"reference_resolved": 0.0,
62
  "build_compatible": 100.0}
app/ngs/stages/stage5_alignment.py CHANGED
@@ -15,9 +15,6 @@ feeds Stage 6 (BAM processing), Stage 7 (alignment QC) and Stage 8 (coverage).
15
 
16
  from __future__ import annotations
17
 
18
- import os
19
- import shutil
20
- import subprocess
21
  from typing import Optional
22
 
23
  from app.ngs.contracts import QcStatus, StageContract, ThresholdRule
@@ -38,10 +35,6 @@ def choose_aligner(assay: str, read_length: Optional[int]) -> str:
38
  return "bwa-mem2"
39
 
40
 
41
- def _find_tool(name: str) -> Optional[str]:
42
- return shutil.which(name)
43
-
44
-
45
  def _align_python(
46
  ref_seq: str,
47
  reads: list[tuple[str, str, str]],
@@ -59,21 +52,15 @@ def align(
59
  ) -> tuple[list[dict], dict]:
60
  """Align reads to a reference sequence. Returns (SAM records, meta).
61
 
62
- External tools are used when available; otherwise the pure-Python fallback runs. The
63
- selection and execution record both go into the provenance/meta so it's auditable.
64
  """
65
  aligner = choose_aligner(assay, read_length)
66
- used_external = False
67
- tool = _find_tool(aligner)
68
- if tool:
69
- used_external = True
70
- # Real BWA/STAR would consume FASTQ + reference FASTA; for the surrogate we still
71
- # produce records from the Python aligner but record the intended external tool.
72
  records = _align_python(ref_seq, reads, ref_name=ref_name)
73
  return records, {
74
- "aligner": aligner,
75
- "aligner_available": bool(tool),
76
- "executor": "external" if used_external else "python-fallback",
77
  "reference": ref_name,
78
  "read_length": read_length,
79
  "aligned_reads": len(records),
@@ -102,7 +89,7 @@ def _stage5_run(sample: dict, state: dict) -> tuple[dict, dict]:
102
  def stage5_contract() -> StageContract:
103
  return StageContract(
104
  step="alignment",
105
- tool="platform-aligner (bwa-mem2 / minimap2 / STAR)",
106
  version="0.1.0",
107
  inputs=["clean_fastq", "reference_sequence"],
108
  outputs=["aligned_reads"],
@@ -112,6 +99,7 @@ def stage5_contract() -> StageContract:
112
  ],
113
  fail_blocks=False,
114
  run=_stage5_run,
 
115
  )
116
 
117
 
 
15
 
16
  from __future__ import annotations
17
 
 
 
 
18
  from typing import Optional
19
 
20
  from app.ngs.contracts import QcStatus, StageContract, ThresholdRule
 
35
  return "bwa-mem2"
36
 
37
 
 
 
 
 
38
  def _align_python(
39
  ref_seq: str,
40
  reads: list[tuple[str, str, str]],
 
52
  ) -> tuple[list[dict], dict]:
53
  """Align reads to a reference sequence. Returns (SAM records, meta).
54
 
55
+ This endpoint currently executes the platform's seed-based Python aligner. The recommended
56
+ production aligner is recorded separately and must never be represented as executed.
57
  """
58
  aligner = choose_aligner(assay, read_length)
 
 
 
 
 
 
59
  records = _align_python(ref_seq, reads, ref_name=ref_name)
60
  return records, {
61
+ "recommended_production_aligner": aligner,
62
+ "executed_implementation": "bionexus-seed-aligner",
63
+ "executor": "python-surrogate",
64
  "reference": ref_name,
65
  "read_length": read_length,
66
  "aligned_reads": len(records),
 
89
  def stage5_contract() -> StageContract:
90
  return StageContract(
91
  step="alignment",
92
+ tool="bionexus-seed-aligner",
93
  version="0.1.0",
94
  inputs=["clean_fastq", "reference_sequence"],
95
  outputs=["aligned_reads"],
 
99
  ],
100
  fail_blocks=False,
101
  run=_stage5_run,
102
+ evidence_level="SURROGATE",
103
  )
104
 
105
 
app/ngs/stages/stage6_bam.py CHANGED
@@ -104,6 +104,7 @@ def stage6_contract() -> StageContract:
104
  ],
105
  fail_blocks=False, # high duplication is a WARN, not a stop (it can be filtered)
106
  run=_stage6_run,
 
107
  )
108
 
109
 
 
104
  ],
105
  fail_blocks=False, # high duplication is a WARN, not a stop (it can be filtered)
106
  run=_stage6_run,
107
+ evidence_level="SURROGATE",
108
  )
109
 
110
 
app/ngs/stages/stage8_coverage.py CHANGED
@@ -156,12 +156,37 @@ def _stage8_run(sample: dict, state: dict) -> tuple[dict, dict]:
156
  eng = coverage_engine(records, ref_lengths, targets)
157
  state.setdefault("coverage", {})["engine"] = eng
158
  g = eng["genome"]
 
 
 
 
 
159
  return eng, {
160
  "coverage_ok": g["coverage_30x"],
161
  "uniformity_ok": g["uniformity"],
162
  }
163
 
164
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
165
  def stage8_contract() -> StageContract:
166
  return StageContract(
167
  step="coverage",
@@ -169,12 +194,7 @@ def stage8_contract() -> StageContract:
169
  version="0.1.0",
170
  inputs=["processed_bam", "reference_lengths", "target_regions?"],
171
  outputs=["coverage_report"],
172
- rules=[
173
- ThresholdRule(name="coverage_ok", metric="coverage_ok",
174
- evaluate=lambda v: _pct_rule(v, 85, 60)),
175
- ThresholdRule(name="uniformity_ok", metric="uniformity_ok",
176
- evaluate=lambda v: _pct_rule(v, 80, 60)),
177
- ],
178
  fail_blocks=False,
179
  run=_stage8_run,
180
  )
 
156
  eng = coverage_engine(records, ref_lengths, targets)
157
  state.setdefault("coverage", {})["engine"] = eng
158
  g = eng["genome"]
159
+ if sample.get("synthetic_reference") or sample.get("demonstration_data"):
160
+ return eng, {
161
+ "synthetic_reference_covered": g["coverage_1x"],
162
+ "synthetic_uniformity": g["uniformity"],
163
+ }
164
  return eng, {
165
  "coverage_ok": g["coverage_30x"],
166
  "uniformity_ok": g["uniformity"],
167
  }
168
 
169
 
170
+ def _coverage_rules(sample: dict) -> list[ThresholdRule]:
171
+ if sample.get("synthetic_reference") or sample.get("demonstration_data"):
172
+ return [
173
+ ThresholdRule(
174
+ name="synthetic_reference_covered", metric="synthetic_reference_covered",
175
+ evaluate=lambda v: _pct_rule(v, 95, 80), expectation=">= 95% at 1x",
176
+ ),
177
+ ThresholdRule(
178
+ name="synthetic_uniformity", metric="synthetic_uniformity",
179
+ evaluate=lambda v: _pct_rule(v, 95, 80), expectation=">= 95%",
180
+ ),
181
+ ]
182
+ return [
183
+ ThresholdRule(name="coverage_ok", metric="coverage_ok",
184
+ evaluate=lambda v: _pct_rule(v, 85, 60), expectation=">= 85% at 30x"),
185
+ ThresholdRule(name="uniformity_ok", metric="uniformity_ok",
186
+ evaluate=lambda v: _pct_rule(v, 80, 60), expectation=">= 80%"),
187
+ ]
188
+
189
+
190
  def stage8_contract() -> StageContract:
191
  return StageContract(
192
  step="coverage",
 
194
  version="0.1.0",
195
  inputs=["processed_bam", "reference_lengths", "target_regions?"],
196
  outputs=["coverage_report"],
197
+ rules=_coverage_rules,
 
 
 
 
 
198
  fail_blocks=False,
199
  run=_stage8_run,
200
  )
app/ngs/stages/stage9_contamination.py CHANGED
@@ -79,14 +79,14 @@ def _stage9_run(sample: dict, state: dict) -> tuple[dict, dict]:
79
  if records is None:
80
  return ({"error": "contamination needs aligned records",
81
  "unevaluated": "no aligned reads"},
82
- {"contam_checked": 0.0, "contam_ok": 100.0})
83
  if not sites:
84
  return ({"unevaluated": "no SNP panel supplied; contamination not asserted",
85
- "contam_rate": 0.0},
86
- {"contam_checked": 0.0, "contam_ok": 100.0})
87
  eng = contamination_engine(records, sites)
88
  state.setdefault("contamination", {})["engine"] = eng
89
- return eng, {"contam_checked": 100.0, "contam_ok": round(100.0 - eng["contam_rate"], 3)}
90
 
91
 
92
  def stage9_contract() -> StageContract:
@@ -97,30 +97,27 @@ def stage9_contract() -> StageContract:
97
  inputs=["processed_bam", "snp_sites"],
98
  outputs=["contamination_report"],
99
  rules=[
100
- ThresholdRule(name="contam_checked", metric="contam_checked",
101
- evaluate=lambda v: _checked_rule(v)),
102
- ThresholdRule(name="contam_ok", metric="contam_ok",
103
- evaluate=lambda v: _invert_pct_rule(v, 98.0, 95.0)),
 
 
104
  ],
105
  fail_blocks=True, # contamination FAIL blocks downstream variant calling
106
  run=_stage9_run,
 
107
  )
108
 
109
 
110
- def _checked_rule(v):
111
- # 100 = a real SNP-panel check ran; 0 = unevaluated (WARN, not blocking)
112
- return QcStatus.PASS if float(v) >= 100.0 else QcStatus.WARN
113
-
114
-
115
- def _invert_pct_rule(v, ok, warn):
116
- # ok = min "purity" for PASS, warn = min for WARN; below warn -> FAIL
117
  try:
118
  v = float(v)
119
  except (TypeError, ValueError):
120
  return QcStatus.FAIL
121
- if v >= ok:
122
  return QcStatus.PASS
123
- if v >= warn:
124
  return QcStatus.WARN
125
  return QcStatus.FAIL
126
 
@@ -131,8 +128,7 @@ def run_contamination(
131
  ) -> dict:
132
  from app.ngs.contracts import apply_rules, QcResult
133
  eng = contamination_engine(records, snp_sites)
134
- metrics = {"contam_ok": round(100.0 - eng["contam_rate"], 3),
135
- "contam_checked": 100.0 if snp_sites else 0.0}
136
  contract = stage9_contract()
137
  result = QcResult.from_metrics(apply_rules(contract.resolve_rules({}), metrics),
138
  fail_blocks=True)
 
79
  if records is None:
80
  return ({"error": "contamination needs aligned records",
81
  "unevaluated": "no aligned reads"},
82
+ {})
83
  if not sites:
84
  return ({"unevaluated": "no SNP panel supplied; contamination not asserted",
85
+ "contam_rate": None, "status": "NOT_EVALUATED"},
86
+ {})
87
  eng = contamination_engine(records, sites)
88
  state.setdefault("contamination", {})["engine"] = eng
89
+ return eng, {"contamination_pct": eng["contam_rate"]}
90
 
91
 
92
  def stage9_contract() -> StageContract:
 
97
  inputs=["processed_bam", "snp_sites"],
98
  outputs=["contamination_report"],
99
  rules=[
100
+ ThresholdRule(
101
+ name="contamination_pct", metric="contamination_pct",
102
+ evaluate=lambda v: _contamination_rule(v), expectation="<= 2%",
103
+ optional=True,
104
+ missing_detail="Not evaluated: a reference-matched SNP panel was not supplied.",
105
+ ),
106
  ],
107
  fail_blocks=True, # contamination FAIL blocks downstream variant calling
108
  run=_stage9_run,
109
+ evidence_level="SURROGATE",
110
  )
111
 
112
 
113
+ def _contamination_rule(v):
 
 
 
 
 
 
114
  try:
115
  v = float(v)
116
  except (TypeError, ValueError):
117
  return QcStatus.FAIL
118
+ if v <= 2.0:
119
  return QcStatus.PASS
120
+ if v <= 5.0:
121
  return QcStatus.WARN
122
  return QcStatus.FAIL
123
 
 
128
  ) -> dict:
129
  from app.ngs.contracts import apply_rules, QcResult
130
  eng = contamination_engine(records, snp_sites)
131
+ metrics = {"contamination_pct": eng["contam_rate"]} if snp_sites else {}
 
132
  contract = stage9_contract()
133
  result = QcResult.from_metrics(apply_rules(contract.resolve_rules({}), metrics),
134
  fail_blocks=True)
app/routers/ngs_v2.py CHANGED
@@ -16,17 +16,31 @@ import logging
16
  import os
17
  import random
18
  import tempfile
 
 
19
  from typing import Optional
20
 
21
- from fastapi import APIRouter, HTTPException
22
  from pydantic import BaseModel
23
 
24
  from app.ngs.assays import AssayRouter
25
  from app.ngs.orchestrator import build_dag, wgs_wes_germline_stages
 
 
 
26
  from app.ngs.visualization import build_visualization
 
 
 
 
 
 
 
 
27
 
28
  logger = logging.getLogger(__name__)
29
  router = APIRouter(prefix="/api/ngs/v2", tags=["ngs-v2"])
 
30
 
31
  _DETECTOR = AssayRouter()
32
 
@@ -101,7 +115,7 @@ class DetectRequest(BaseModel):
101
  metadata: dict = {}
102
 
103
 
104
- def _read_fastq(path: str, cap: int = 2000) -> list[tuple[str, str, str]]:
105
  if not os.path.isfile(path):
106
  raise HTTPException(status_code=400, detail=f"file not found: {path}")
107
  reads: list[tuple[str, str, str]] = []
@@ -118,9 +132,10 @@ def _read_fastq(path: str, cap: int = 2000) -> list[tuple[str, str, str]]:
118
  if not name.startswith("@") or not plus.startswith("+"):
119
  continue
120
  reads.append((name[1:].split()[0], seq, qual))
 
121
  except Exception as exc:
122
  raise HTTPException(status_code=400, detail=f"failed to read {path}: {exc}")
123
- return reads
124
 
125
 
126
  def _demo_sequence(rng: random.Random, length: int = 150) -> str:
@@ -196,7 +211,8 @@ def _detection(payload, files: Optional[list[str]] = None) -> dict:
196
  def list_stages():
197
  contracts = wgs_wes_germline_stages()
198
  return {"pipeline": "wgs-wes-germline", "stages": [
199
- {"step": c.step, "tool": c.tool, "inputs": c.inputs, "outputs": c.outputs,
 
200
  "fail_blocks": c.fail_blocks, "expectation": _STAGE_INTRO.get(c.step, "")}
201
  for c in contracts
202
  ]}
@@ -207,6 +223,66 @@ def list_demos():
207
  return {"demos": [{"id": key, **value} for key, value in _DEMO_PROFILES.items()]}
208
 
209
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
210
  @router.post("/detect")
211
  def detect(payload: DetectRequest):
212
  return _detection(payload)
@@ -222,9 +298,12 @@ def analyze(payload: AnalyzeRequest):
222
  raise HTTPException(status_code=400, detail="provide file_paths, fastq_url, or demo_profile")
223
 
224
  reads_all: dict[str, list[tuple[str, str, str]]] = {}
 
225
  for f in files:
226
- got = _read_fastq(f)
227
  reads_all[f] = got
 
 
228
  if not got:
229
  raise HTTPException(status_code=400, detail=f"no FASTQ records read from {f}")
230
 
@@ -233,6 +312,8 @@ def analyze(payload: AnalyzeRequest):
233
  assay = (payload.assay or demo_assay or detection["assay"] or "WGS").upper()
234
  if assay in ("UNKNOWN", ""):
235
  assay = "WGS"
 
 
236
 
237
  combined_reads = [read for f in files for read in reads_all.get(f, [])]
238
  sample: dict = {
@@ -240,7 +321,11 @@ def analyze(payload: AnalyzeRequest):
240
  "reference": payload.reference or "grch38",
241
  "assay": assay,
242
  "sample_type": payload.sample_type or detection["sample_type"],
243
- "metadata": {**(payload.metadata or {}), **({"demo_profile": payload.demo_profile} if payload.demo_profile else {})},
 
 
 
 
244
  "demonstration_data": bool(payload.demo_profile or (payload.metadata or {}).get("demonstration_data")),
245
  "synthetic_reference": bool(payload.synthetic_reference or payload.demo_profile),
246
  "reads": combined_reads,
@@ -252,15 +337,30 @@ def analyze(payload: AnalyzeRequest):
252
 
253
  pipe = build_dag(assay)
254
  report = pipe.run(sample)
 
 
 
 
 
 
 
 
 
 
 
255
  return {
256
  "detection": detection,
257
  "requested": {
258
  "assay": assay,
259
- "reference": sample["reference"],
 
260
  "synthetic_reference": use_synthetic_reference,
261
  "demo_profile": payload.demo_profile,
262
  "reads_loaded": {os.path.basename(f): len(reads_all[f]) for f in files},
263
  "reads_analyzed": len(combined_reads),
 
 
 
264
  },
265
  "demo": demo,
266
  "pipeline": report,
 
16
  import os
17
  import random
18
  import tempfile
19
+ import json
20
+ from pathlib import Path
21
  from typing import Optional
22
 
23
+ from fastapi import APIRouter, Depends, HTTPException
24
  from pydantic import BaseModel
25
 
26
  from app.ngs.assays import AssayRouter
27
  from app.ngs.orchestrator import build_dag, wgs_wes_germline_stages
28
+ from app.ngs.production import build_production_plan, evaluate_clinical_evidence
29
+ from app.ngs.execution import artifact_manifest, executor_capabilities, get_run, submit_run
30
+ from app.services.auth import require_user_id
31
  from app.ngs.visualization import build_visualization
32
+ from app.models.responses import (
33
+ NgsClinicalEvidenceRequest,
34
+ NgsClinicalEvidenceResponse,
35
+ NgsProductionPlanRequest,
36
+ NgsProductionPlanResponse,
37
+ NgsProductionRunResponse,
38
+ NgsProductionSubmitResponse,
39
+ )
40
 
41
  logger = logging.getLogger(__name__)
42
  router = APIRouter(prefix="/api/ngs/v2", tags=["ngs-v2"])
43
+ FASTQ_RECORD_CAP_PER_FILE = 2000
44
 
45
  _DETECTOR = AssayRouter()
46
 
 
115
  metadata: dict = {}
116
 
117
 
118
+ def _read_fastq(path: str, cap: int = FASTQ_RECORD_CAP_PER_FILE) -> tuple[list[tuple[str, str, str]], bool]:
119
  if not os.path.isfile(path):
120
  raise HTTPException(status_code=400, detail=f"file not found: {path}")
121
  reads: list[tuple[str, str, str]] = []
 
132
  if not name.startswith("@") or not plus.startswith("+"):
133
  continue
134
  reads.append((name[1:].split()[0], seq, qual))
135
+ truncated = bool(fh.readline())
136
  except Exception as exc:
137
  raise HTTPException(status_code=400, detail=f"failed to read {path}: {exc}")
138
+ return reads, truncated
139
 
140
 
141
  def _demo_sequence(rng: random.Random, length: int = 150) -> str:
 
211
  def list_stages():
212
  contracts = wgs_wes_germline_stages()
213
  return {"pipeline": "wgs-wes-germline", "stages": [
214
+ {"step": c.step, "operation": c.step.replace("_", " "), "implementation": c.tool,
215
+ "evidence_level": c.evidence_level, "inputs": c.inputs, "outputs": c.outputs,
216
  "fail_blocks": c.fail_blocks, "expectation": _STAGE_INTRO.get(c.step, "")}
217
  for c in contracts
218
  ]}
 
223
  return {"demos": [{"id": key, **value} for key, value in _DEMO_PROFILES.items()]}
224
 
225
 
226
+ @router.get("/benchmarks/portable")
227
+ def portable_benchmark():
228
+ report_path = Path(__file__).resolve().parents[1] / "ngs" / "portable_benchmark_report.json"
229
+ with report_path.open("r", encoding="utf-8") as handle:
230
+ return json.load(handle)
231
+
232
+
233
+ @router.post("/production/plan", response_model=NgsProductionPlanResponse)
234
+ def production_plan(payload: NgsProductionPlanRequest):
235
+ return build_production_plan(payload)
236
+
237
+
238
+ @router.get("/production/capabilities")
239
+ def production_capabilities():
240
+ return executor_capabilities()
241
+
242
+
243
+ @router.post("/production/submit", response_model=NgsProductionSubmitResponse)
244
+ def production_submit(payload: NgsProductionPlanRequest, user_id: str = Depends(require_user_id)):
245
+ plan = build_production_plan(payload)
246
+ if not plan["ready_to_launch"]:
247
+ raise HTTPException(status_code=422, detail={"message": "production launch contract is blocked", "blockers": plan["blockers"]})
248
+ executor = "awsbatch" if payload.execution_profile == "awsbatch" else "slurm" if payload.execution_profile == "slurm" else "local"
249
+ try:
250
+ run = submit_run(executor, plan["command_argv"], payload.outdir, user_id)
251
+ except RuntimeError as exc:
252
+ raise HTTPException(status_code=503, detail=str(exc)) from exc
253
+ except Exception as exc:
254
+ logger.exception("Production Sarek submission failed")
255
+ raise HTTPException(status_code=502, detail=f"executor submission failed: {type(exc).__name__}") from exc
256
+ return {
257
+ "run_id": run["run_id"], "state": "SUBMITTED", "executor": executor,
258
+ "executor_job_id": run["executor_job_id"], "message": "Real nf-core/sarek run submitted; no exploratory fallback was used.",
259
+ }
260
+
261
+
262
+ @router.get("/production/runs/{run_id}", response_model=NgsProductionRunResponse)
263
+ def production_run_status(run_id: str, user_id: str = Depends(require_user_id)):
264
+ try:
265
+ return get_run(run_id, user_id)
266
+ except FileNotFoundError as exc:
267
+ raise HTTPException(status_code=404, detail="production run not found") from exc
268
+
269
+
270
+ @router.get("/production/runs/{run_id}/artifacts")
271
+ def production_run_artifacts(run_id: str, user_id: str = Depends(require_user_id)):
272
+ try:
273
+ return artifact_manifest(run_id, user_id)
274
+ except FileNotFoundError as exc:
275
+ raise HTTPException(status_code=404, detail="production run not found") from exc
276
+ except Exception as exc:
277
+ logger.exception("Production Sarek artifact import failed")
278
+ raise HTTPException(status_code=502, detail=f"artifact import failed: {type(exc).__name__}") from exc
279
+
280
+
281
+ @router.post("/clinical/evaluate", response_model=NgsClinicalEvidenceResponse)
282
+ def clinical_evidence(payload: NgsClinicalEvidenceRequest):
283
+ return evaluate_clinical_evidence(payload)
284
+
285
+
286
  @router.post("/detect")
287
  def detect(payload: DetectRequest):
288
  return _detection(payload)
 
298
  raise HTTPException(status_code=400, detail="provide file_paths, fastq_url, or demo_profile")
299
 
300
  reads_all: dict[str, list[tuple[str, str, str]]] = {}
301
+ truncated_inputs: list[str] = []
302
  for f in files:
303
+ got, truncated = _read_fastq(f)
304
  reads_all[f] = got
305
+ if truncated:
306
+ truncated_inputs.append(os.path.basename(f))
307
  if not got:
308
  raise HTTPException(status_code=400, detail=f"no FASTQ records read from {f}")
309
 
 
312
  assay = (payload.assay or demo_assay or detection["assay"] or "WGS").upper()
313
  if assay in ("UNKNOWN", ""):
314
  assay = "WGS"
315
+ if payload.demo_profile:
316
+ detection["sample_type"] = "synthetic-positive-control"
317
 
318
  combined_reads = [read for f in files for read in reads_all.get(f, [])]
319
  sample: dict = {
 
321
  "reference": payload.reference or "grch38",
322
  "assay": assay,
323
  "sample_type": payload.sample_type or detection["sample_type"],
324
+ "metadata": {
325
+ **({"platform": "illumina-simulated"} if payload.demo_profile else {}),
326
+ **(payload.metadata or {}),
327
+ **({"demo_profile": payload.demo_profile} if payload.demo_profile else {}),
328
+ },
329
  "demonstration_data": bool(payload.demo_profile or (payload.metadata or {}).get("demonstration_data")),
330
  "synthetic_reference": bool(payload.synthetic_reference or payload.demo_profile),
331
  "reads": combined_reads,
 
337
 
338
  pipe = build_dag(assay)
339
  report = pipe.run(sample)
340
+ report["validation"]["input_sampling"] = {
341
+ "mode": "SAMPLED_PREVIEW" if truncated_inputs else "COMPLETE_FOR_SUPPLIED_FILES",
342
+ "record_cap_per_file": FASTQ_RECORD_CAP_PER_FILE,
343
+ "truncated_files": truncated_inputs,
344
+ "all_records_processed": not truncated_inputs,
345
+ }
346
+ if truncated_inputs:
347
+ report["warnings"].append(
348
+ f"Only the first {FASTQ_RECORD_CAP_PER_FILE} records per FASTQ were analyzed; "
349
+ "this result is a sampled preview and cannot support whole-run QC or variant conclusions."
350
+ )
351
  return {
352
  "detection": detection,
353
  "requested": {
354
  "assay": assay,
355
+ "reference": "synthetic-positive-control" if use_synthetic_reference else sample["reference"],
356
+ "reference_template_requested": sample["reference"] if use_synthetic_reference else None,
357
  "synthetic_reference": use_synthetic_reference,
358
  "demo_profile": payload.demo_profile,
359
  "reads_loaded": {os.path.basename(f): len(reads_all[f]) for f in files},
360
  "reads_analyzed": len(combined_reads),
361
+ "record_cap_per_file": FASTQ_RECORD_CAP_PER_FILE,
362
+ "truncated_files": truncated_inputs,
363
+ "all_records_processed": not truncated_inputs,
364
  },
365
  "demo": demo,
366
  "pipeline": report,
app/routers/pipeline_v2.py CHANGED
@@ -1,1108 +1,1108 @@
1
- """
2
- In-memory pipeline v2 — runs BLAST → UniProt → MSA → Phylo → Domains → Interpretation
3
- in a background thread. Uses a thread-safe dict for job storage.
4
- """
5
-
6
- import asyncio
7
- import logging
8
- import threading
9
- import uuid
10
- from datetime import datetime, timezone
11
-
12
- from fastapi import APIRouter, HTTPException, Depends, Request
13
- from pydantic import BaseModel, Field
14
-
15
- from app.config import settings
16
- from app.deps import limiter
17
- from app.services.rate_limit import check_daily_limit_pipelines
18
- from app.services.auth import get_user_id
19
- from app.services.supabase import get_supabase
20
- from app.integrations.ncbi import blast as ncbi_blast
21
- from app.integrations.ncbi.parser import parse_blast_xml
22
- from app.services.validators import validate_fasta
23
- from app.services.sequence_utils import detect_source_from_accession, map_refseq_to_uniprot, detect_sequence_type
24
- from app.services.blast_config import resolve_blast_params
25
- from app.services.identifier_resolution import resolve_to_uniprot
26
- from app.tools.ebi_msa import EBI_TOOLS, run_ebi_msa, run_ebi_msa_best_effort
27
- from app.tools.blast import BlastTool
28
- from app.tools.uniprot import UniprotTool
29
- from app.tools.pairwise_alignment import pairwise_align, VALID_MODES
30
-
31
- logger = logging.getLogger(__name__)
32
- router = APIRouter()
33
-
34
- _jobs: dict[str, dict] = {}
35
- _jobs_lock = threading.Lock()
36
-
37
- STEP_ORDER = ["blast", "uniprot", "msa", "phylo", "domains", "pathway_enrichment", "alphafold", "interpret"]
38
-
39
-
40
- def _get_job(job_id: str) -> dict | None:
41
- with _jobs_lock:
42
- return _jobs.get(job_id)
43
-
44
-
45
- def _is_real_uuid(job_id: str) -> bool:
46
- try:
47
- uuid.UUID(job_id)
48
- return True
49
- except (ValueError, AttributeError):
50
- return False
51
-
52
-
53
- def _persist_v2_final(job_id: str, status: str, context: dict, error: str | None = None):
54
- """Update the persisted wizard job row with the final status + full context."""
55
- if not _is_real_uuid(job_id):
56
- return
57
- now = datetime.now(timezone.utc).isoformat()
58
- payload = {
59
- "status": status,
60
- "progress_pct": 100,
61
- "completed_at": now,
62
- "context_json": context,
63
- "error_message": error,
64
- }
65
- if error:
66
- payload["error"] = error
67
- try:
68
- get_supabase().table("jobs").update(payload).eq("id", job_id).execute()
69
- except Exception as e:
70
- logger.warning("Could not update v2 job %s to Supabase: %s", job_id, e)
71
-
72
-
73
- def _persist_v2_job(job_id: str, payload: dict):
74
- """Persist a wizard v2 job to Supabase (best-effort) so it shows up in
75
- the user's job history and can be shared. Uses the service role key."""
76
- try:
77
- get_supabase().table("jobs").insert(payload).execute()
78
- except Exception as e:
79
- logger.warning("Could not persist v2 job %s to Supabase: %s", job_id, e)
80
-
81
-
82
- def _set_step_status(job_id: str, step: str, status: str, progress: int = 0, data: dict | None = None, error: str | None = None):
83
- with _jobs_lock:
84
- if job_id not in _jobs:
85
- return
86
- _jobs[job_id]["steps"][step] = {"status": status, "progress": progress, "data": data, "error": error}
87
- if status == "running":
88
- _jobs[job_id]["current_step"] = step
89
-
90
-
91
- def _set_job_failed(job_id: str, message: str):
92
- with _jobs_lock:
93
- if job_id in _jobs:
94
- _jobs[job_id]["status"] = "failed"
95
- _jobs[job_id]["error"] = message
96
-
97
-
98
- class PipelineV2RunRequest(BaseModel):
99
- sequence: str = Field(..., min_length=6, description="Protein sequence (FASTA or raw)")
100
- steps: list[str] = Field(default_factory=lambda: list(STEP_ORDER), description="Steps to run")
101
- fast_mode: bool = Field(default=False, description="Use Swiss-Prot instead of nr for faster results")
102
- database: str = Field("", description="BLAST database override")
103
- program: str = Field("", description="BLAST program override")
104
- max_hits: int = Field(100, description="Max BLAST hits to return")
105
- query_accession: str = Field("", description="Optional query accession for display")
106
- alignment_mode: str = Field("global", description="Alignment mode for the MSA step: 'global' (full-length) or 'local' (Smith-Waterman refinement of query vs top hit)")
107
- parent_job_id: str | None = Field(None, description="Parent job ID for DAG branching")
108
-
109
-
110
- @router.post("/run")
111
- async def run_pipeline_v2(request: Request, req: PipelineV2RunRequest, user_id: str | None = Depends(get_user_id)):
112
- validation = validate_fasta(req.sequence, "blast")
113
- if not validation.valid:
114
- raise HTTPException(status_code=400, detail=validation.error)
115
-
116
- alignment_mode = req.alignment_mode.lower()
117
- if alignment_mode not in VALID_MODES:
118
- raise HTTPException(status_code=400, detail=f"alignment_mode must be one of {list(VALID_MODES)}, got {alignment_mode!r}")
119
-
120
- seq = str(validation.sequences[0].seq).upper()
121
- clean = "".join(c for c in seq if c.isalpha())
122
-
123
- job_id = str(uuid.uuid4())
124
- now = datetime.now(timezone.utc).isoformat()
125
-
126
- requested = [s for s in req.steps if s in STEP_ORDER]
127
- if not requested:
128
- requested = list(STEP_ORDER)
129
-
130
- steps_dict = {s: {"status": "pending", "progress": 0, "data": None, "error": None} for s in STEP_ORDER}
131
-
132
- blast_params = {
133
- "database": req.database,
134
- "program": req.program,
135
- "max_hits": req.max_hits,
136
- "query_accession": req.query_accession,
137
- }
138
-
139
- with _jobs_lock:
140
- _jobs[job_id] = {
141
- "job_id": job_id,
142
- "status": "running",
143
- "current_step": None,
144
- "steps": steps_dict,
145
- "requested_steps": requested,
146
- "sequence": clean,
147
- "blast_params": blast_params,
148
- "alignment_mode": alignment_mode,
149
- "error": None,
150
- "created_at": now,
151
- }
152
-
153
- # Best-effort persistence so wizard jobs appear in history and can be shared.
154
- persist_payload = {
155
- "id": job_id,
156
- "user_id": user_id,
157
- "tool": "wizard_v2",
158
- "pipeline_type": "wizard_v2",
159
- "query_preview": clean[:60],
160
- "status": "running",
161
- "progress_pct": 0,
162
- "title": "Wizard pipeline",
163
- "description": "Pipeline v2 wizard run",
164
- "created_at": now,
165
- }
166
- if req.parent_job_id:
167
- persist_payload["parent_job_id"] = req.parent_job_id
168
- _persist_v2_job(job_id, persist_payload)
169
-
170
- t = threading.Thread(
171
- target=_run_pipeline,
172
- args=(job_id, clean, requested),
173
- kwargs={"fast_mode": req.fast_mode, "blast_params": blast_params, "alignment_mode": alignment_mode},
174
- daemon=True,
175
- )
176
- t.start()
177
-
178
- return {"job_id": job_id}
179
-
180
-
181
- @router.get("/status/{job_id}")
182
- async def get_pipeline_v2_status(job_id: str):
183
- job = _get_job(job_id)
184
- if not job:
185
- raise HTTPException(status_code=404, detail="Job not found")
186
- return job
187
-
188
-
189
- async def run_pipeline(
190
- sequence: str,
191
- organism: str = "Homo sapiens",
192
- analysis_type: str = "comprehensive",
193
- status_callback=None,
194
- fast_mode: bool = False,
195
- blast_params: dict | None = None,
196
- ) -> dict:
197
- """Public async entry point for the pipeline (used by pipeline_worker).
198
-
199
- Creates a temporary in-memory job, runs the configured steps, and
200
- returns the context dict with all results.
201
- """
202
- job_id = f"worker-{uuid.uuid4().hex[:12]}"
203
- requested = list(STEP_ORDER)
204
- steps_dict = {s: {"status": "pending", "progress": 0, "data": None, "error": None} for s in STEP_ORDER}
205
-
206
- with _jobs_lock:
207
- _jobs[job_id] = {
208
- "job_id": job_id,
209
- "status": "running",
210
- "current_step": None,
211
- "steps": steps_dict,
212
- "requested_steps": requested,
213
- "sequence": sequence,
214
- "blast_params": blast_params or {},
215
- "error": None,
216
- "created_at": datetime.now(timezone.utc).isoformat(),
217
- }
218
-
219
- try:
220
- await _execute(
221
- job_id,
222
- sequence,
223
- requested,
224
- status_callback=status_callback,
225
- fast_mode=fast_mode,
226
- blast_params=blast_params,
227
- )
228
- finally:
229
- job = _get_job(job_id)
230
- with _jobs_lock:
231
- _jobs.pop(job_id, None)
232
-
233
- if job and job.get("status") == "failed":
234
- raise RuntimeError(job.get("error", "Pipeline failed"))
235
-
236
- query_accession = ((blast_params or {}).get("query_accession") or "").strip()
237
- context: dict = {
238
- "sequence": sequence,
239
- "length": len(sequence),
240
- "query": {
241
- "sequence": sequence,
242
- "length": len(sequence),
243
- "sequence_type": detect_sequence_type(sequence) or "protein",
244
- },
245
- }
246
- if query_accession:
247
- context["query"]["accession"] = query_accession
248
- if job:
249
- for step_name, step_info in job.get("steps", {}).items():
250
- if step_info.get("data"):
251
- context[step_name] = step_info["data"]
252
- return context
253
-
254
-
255
- # ---------------------------------------------------------------------------
256
- # Background pipeline
257
- # ---------------------------------------------------------------------------
258
-
259
- def _run_pipeline(job_id: str, sequence: str, steps: list[str], fast_mode: bool = False, blast_params: dict | None = None, alignment_mode: str = "global"):
260
- loop = asyncio.new_event_loop()
261
- asyncio.set_event_loop(loop)
262
- try:
263
- loop.run_until_complete(
264
- _execute(job_id, sequence, steps, fast_mode=fast_mode, blast_params=blast_params, alignment_mode=alignment_mode)
265
- )
266
- except Exception as e:
267
- logger.exception(f"[{job_id}] Unhandled pipeline error")
268
- _set_job_failed(job_id, f"Pipeline error: {e}")
269
- finally:
270
- loop.close()
271
- asyncio.set_event_loop(None)
272
-
273
-
274
- async def _execute(job_id: str, sequence: str, steps: list[str], status_callback=None, fast_mode: bool = False, blast_params: dict | None = None, alignment_mode: str = "global"):
275
- context: dict = {
276
- "sequence": sequence,
277
- "length": len(sequence),
278
- "query": {
279
- "sequence": sequence,
280
- "length": len(sequence),
281
- "sequence_type": detect_sequence_type(sequence) or "protein",
282
- },
283
- }
284
-
285
- _STEP_FRONTEND = {
286
- "blast": "running",
287
- "uniprot": "fetching_uniprot",
288
- "msa": "running_msa",
289
- "phylo": "running_msa",
290
- "domains": "fetching_uniprot",
291
- "pathway_enrichment": "pathway_enrichment",
292
- "alphafold": "fetching_alphafold",
293
- "interpret": "interpreting",
294
- }
295
-
296
- _failed_step = None
297
- _failed_error = None
298
-
299
- async def _notify(step_key: str):
300
- if status_callback:
301
- try:
302
- await status_callback(_STEP_FRONTEND.get(step_key, "running"))
303
- except Exception:
304
- pass
305
-
306
- def _mark(step_key: str, status: str, **kw):
307
- _set_step_status(job_id, step_key, status, **kw)
308
-
309
- def _fail(step_key: str, msg: str):
310
- nonlocal _failed_step, _failed_error
311
- _mark(step_key, "failed", error=msg)
312
- _failed_step = step_key
313
- _failed_error = msg
314
-
315
- # ---- Step 1: BLAST (must run first) ----
316
- denovo_mode = False
317
- if "blast" in steps:
318
- await _notify("blast")
319
- _mark("blast", "running", progress=10)
320
- result = await _run_blast(
321
- sequence,
322
- status_callback=status_callback,
323
- fast_mode=fast_mode,
324
- blast_params=blast_params,
325
- )
326
- zero_hits = result.get("count", 0) == 0
327
- if zero_hits and (detect_sequence_type(sequence) or "protein") == "protein":
328
- # Tier 6: no database match at all — characterize from sequence
329
- # alone instead of failing the run (techspec.md §1).
330
- denovo_mode = True
331
- result["_note"] = "No BLAST hits found — switching to de novo characterization"
332
- _mark("blast", "complete", progress=100, data=result)
333
- else:
334
- _mark("blast", "failed" if zero_hits else "complete", progress=100, data=result)
335
- if zero_hits:
336
- _failed_step = "blast"
337
- _failed_error = result.get("error", "No BLAST hits found")
338
- context["blast"] = result
339
-
340
- if denovo_mode:
341
- context["query"]["confidence"] = "de_novo"
342
- await _run_denovo_steps(job_id, sequence, steps, _mark, context)
343
- await _finalize_context(job_id, context)
344
- with _jobs_lock:
345
- if job_id in _jobs:
346
- _jobs[job_id]["status"] = "complete"
347
- _jobs[job_id]["context"] = context
348
- _persist_v2_final(job_id, "complete", context)
349
- return
350
-
351
- # ---- Step 2: Fan-out — UniProt, MSA, Pathway run in parallel ----
352
- # They all only depend on BLAST results, not on each other.
353
- blast_data = context.get("blast", {})
354
- hits = (blast_data.get("hits") if isinstance(blast_data, dict) else []) or []
355
- top_hit = blast_data.get("top_hit") if isinstance(blast_data, dict) else None
356
-
357
- async def _do_uniprot():
358
- candidates = ([top_hit] + hits[:5]) if top_hit else hits[:5]
359
- last_result: dict | None = None
360
- for idx, candidate in enumerate(candidates):
361
- # Only the top candidate may trigger the expensive EBI sequence
362
- # BLAST fallback — otherwise a pathological query could run it up
363
- # to 6 times (up to ~3 min each).
364
- last_result = await _run_uniprot(
365
- candidate,
366
- query_sequence=sequence,
367
- try_sequence=idx == 0,
368
- )
369
- if "error" not in last_result:
370
- return last_result
371
- return last_result if last_result else {"error": "No BLAST hits for UniProt lookup"}
372
-
373
- async def _do_msa():
374
- if not hits:
375
- return {"error": "No BLAST hits for MSA"}
376
- return await _run_msa(sequence, hits, alignment_mode)
377
-
378
- async def _do_pathway():
379
- return await _run_pathway_enrichment(context)
380
-
381
- fan_out = []
382
- fan_names = []
383
- if "uniprot" in steps and not _failed_step:
384
- fan_out.append(_do_uniprot())
385
- fan_names.append("uniprot")
386
- if "msa" in steps and not _failed_step:
387
- fan_out.append(_do_msa())
388
- fan_names.append("msa")
389
- if "pathway_enrichment" in steps and not _failed_step:
390
- fan_out.append(_do_pathway())
391
- fan_names.append("pathway_enrichment")
392
-
393
- if fan_out:
394
- # Notify for the first active step in the fan-out
395
- await _notify(fan_names[0])
396
- for name in fan_names:
397
- _mark(name, "running", progress=10)
398
-
399
- results = await asyncio.gather(*fan_out, return_exceptions=True)
400
-
401
- for name, res in zip(fan_names, results):
402
- if isinstance(res, Exception):
403
- _fail(name, str(res)[:500])
404
- continue
405
-
406
- if name == "uniprot":
407
- s = "complete" if "error" not in res else "failed"
408
- _mark("uniprot", s, progress=100, data=res)
409
- context["uniprot"] = res
410
- if res.get("confidence"):
411
- context["query"]["confidence"] = res["confidence"]
412
- if "error" in res:
413
- _failed_step = "uniprot"
414
- _failed_error = res["error"]
415
-
416
- elif name == "msa":
417
- s = "complete" if res.get("aln_fasta") else "failed"
418
- _mark("msa", s, progress=100, data=res)
419
- context["msa"] = res
420
- if res.get("phylotree"):
421
- context.setdefault("phylo_data", {})["phylotree_newick"] = res["phylotree"]
422
-
423
- elif name == "pathway_enrichment":
424
- s = "complete" if res and res.get("pathways") else "failed"
425
- _mark("pathway_enrichment", s, progress=100, data=res or {})
426
- context["pathway_enrichment"] = res
427
-
428
- # ---- Step 3: Phylo (instant — copies from MSA) ----
429
- if "phylo" in steps and not _failed_step:
430
- _mark("phylo", "running", progress=10)
431
- newick = None
432
- msa_data = context.get("msa", {})
433
- if isinstance(msa_data, dict):
434
- newick = msa_data.get("phylotree")
435
- if not newick:
436
- newick = context.get("phylo_data", {}).get("phylotree_newick")
437
- if newick:
438
- _mark("phylo", "complete", progress=100, data={"phylotree_newick": newick})
439
- context["phylo"] = {"phylotree_newick": newick}
440
- else:
441
- _mark("phylo", "failed", error="No phylotree available from MSA")
442
-
443
- # ---- Step 4: Domains + AlphaFold in parallel ----
444
- # When UniProt resolution exhausted all tiers, these fall back to their
445
- # sequence-only equivalents instead of being silently skipped (§1.2).
446
- uniprot_data = context.get("uniprot", {})
447
- accession = uniprot_data.get("accession") if isinstance(uniprot_data, dict) else None
448
- resolved_uniprot = bool(uniprot_data.get("resolved_uniprot")) if isinstance(uniprot_data, dict) else False
449
-
450
- post_uniprot = []
451
- post_uniprot_names = []
452
- if "domains" in steps and not _failed_step:
453
- post_uniprot.append(
454
- _run_domains_or_denovo(sequence, accession, resolved_uniprot)
455
- )
456
- post_uniprot_names.append("domains")
457
- if "alphafold" in steps and not _failed_step:
458
- post_uniprot.append(
459
- _run_alphafold_or_esmfold(context, sequence, accession, resolved_uniprot)
460
- )
461
- post_uniprot_names.append("alphafold")
462
-
463
- if post_uniprot:
464
- await _notify(post_uniprot_names[0])
465
- for name in post_uniprot_names:
466
- _mark(name, "running", progress=10)
467
-
468
- results2 = await asyncio.gather(*post_uniprot, return_exceptions=True)
469
-
470
- for name, res in zip(post_uniprot_names, results2):
471
- if isinstance(res, Exception):
472
- _fail(name, str(res)[:500])
473
- continue
474
-
475
- if name == "domains":
476
- s = "complete" if res.get("domains") is not None else "failed"
477
- _mark("domains", s, progress=100, data=res)
478
- context["domains"] = res
479
- elif name == "alphafold":
480
- s = "complete" if res else "failed"
481
- _mark("alphafold", s, progress=100, data=res or {})
482
- context["alphafold"] = res
483
-
484
- # ---- Step 5: Interpret (needs all context) ----
485
- if "interpret" in steps and not _failed_step:
486
- await _notify("interpret")
487
- _mark("interpret", "running", progress=10)
488
- result = await _run_interpret(context)
489
- s = "complete" if result.get("interpretation") else "failed"
490
- _mark("interpret", s, progress=100, data=result)
491
- context["interpret"] = result
492
-
493
- # ---- Final status ----
494
- if _failed_step and _failed_step in ("blast", "uniprot"):
495
- with _jobs_lock:
496
- if job_id in _jobs:
497
- _jobs[job_id]["status"] = "failed"
498
- _jobs[job_id]["error"] = f"Pipeline failed at {_failed_step}: {_failed_error}"
499
- _persist_v2_final(job_id, "failed", context, error=f"Pipeline failed at {_failed_step}: {_failed_error}")
500
- else:
501
- await _finalize_context(job_id, context)
502
- with _jobs_lock:
503
- if job_id in _jobs:
504
- _jobs[job_id]["status"] = "complete"
505
- _jobs[job_id]["context"] = context
506
- _persist_v2_final(job_id, "complete", context)
507
-
508
-
509
- async def _finalize_context(job_id: str, context: dict):
510
- """Pre-persist hook: capture source pages (§3) + build the final report."""
511
- user_id = None
512
- with _jobs_lock:
513
- entry = _jobs.get(job_id)
514
- if entry:
515
- user_id = entry.get("user_id")
516
- try:
517
- _capture_run_sources(job_id, context, user_id)
518
- except Exception as e: # captures must never break a run
519
- logger.warning("[%s] page capture wiring failed: %s", job_id, e)
520
- try:
521
- from app.services.final_synthesis import synthesize
522
-
523
- context["final_report"] = await synthesize(context)
524
- except Exception as e:
525
- logger.warning("[%s] final synthesis failed: %s", job_id, e)
526
-
527
-
528
- def _capture_run_sources(job_id: str, context: dict, user_id: str | None):
529
- """Queue one page_captures row per external source this run queried."""
530
- from app.services.page_capture import capture_bg
531
-
532
- blast = context.get("blast") or {}
533
- top_hit = blast.get("top_hit")
534
- if top_hit and top_hit.get("accession"):
535
- capture_bg(job_id, "ncbi", f"https://www.ncbi.nlm.nih.gov/protein/{top_hit['accession']}", user_id)
536
-
537
- uniprot = context.get("uniprot") or {}
538
- if uniprot.get("_de_novo"):
539
- pass # no external annotation source was queried in de novo mode
540
- elif uniprot.get("accession"):
541
- capture_bg(job_id, "uniprot", f"https://www.uniprot.org/uniprotkb/{uniprot['accession']}", user_id)
542
- pdb_ids = uniprot.get("pdb_ids") or []
543
- if pdb_ids:
544
- capture_bg(job_id, "rcsb", f"https://www.rcsb.org/structure/{pdb_ids[0]}", user_id)
545
-
546
- domains = context.get("domains") or {}
547
- dom_list = domains.get("domains") or []
548
- if dom_list and dom_list[0].get("accession"):
549
- acc = dom_list[0]["accession"]
550
- if str(acc).startswith("IPR"):
551
- capture_bg(job_id, "interpro", f"https://www.ebi.ac.uk/interpro/entry/InterPro/{acc}", user_id)
552
-
553
- af = context.get("alphafold") or {}
554
- if af.get("structure_available") and af.get("source") != "esmfold" and uniprot.get("accession"):
555
- capture_bg(job_id, "alphafold", f"https://alphafold.ebi.ac.uk/uniprot/{uniprot['accession']}", user_id)
556
-
557
- pathway = context.get("pathway_enrichment") or {}
558
- pw_list = (pathway.get("pathways") if isinstance(pathway, dict) else None) or []
559
- if pw_list and pw_list[0].get("stId"):
560
- # run_enrichment queries the Reactome projection API exclusively
561
- capture_bg(
562
- job_id, "reactome",
563
- f"https://reactome.org/PathwayBrowser/#{pw_list[0]['stId']}",
564
- user_id,
565
- )
566
-
567
-
568
- # ---------------------------------------------------------------------------
569
- # Step implementations
570
- # ---------------------------------------------------------------------------
571
-
572
- # NCBI database name → EBI ncbiblast database name. EBI has no direct "nr"
573
- # equivalent; uniprotkb is its closest comprehensive protein database.
574
- EBI_BLAST_DATABASE_MAP = {
575
- "nr": "uniprotkb",
576
- "swissprot": "uniprotkb_swissprot",
577
- "pdb": "pdb",
578
- "pdbaa": "pdb_nr",
579
- "refseq_protein": "refseq_protein",
580
- "env_nr": "env_nr",
581
- "nt": "nt",
582
- "refseq_rna": "refseq_rna",
583
- "refseq_genomic": "refseq_genomic",
584
- "est": "est",
585
- "gss": "gss",
586
- }
587
-
588
-
589
- def _build_blast_result(
590
- hits: list[dict],
591
- *,
592
- source: str,
593
- database: str,
594
- program: str,
595
- seq_type: str,
596
- query_accession: str,
597
- query_length: int,
598
- display_limit: int = 20,
599
- ) -> dict:
600
- """Normalize parsed hits (NCBI XML parser or EBI tool) into the pipeline's
601
- canonical BLAST result shape consumed by the frontend."""
602
- top_hit = hits[0] if hits else None
603
- return {
604
- "count": len(hits),
605
- "source": source,
606
- "database": database,
607
- "program": program,
608
- "query_sequence_type": seq_type,
609
- "query_accession": query_accession,
610
- "query_length": query_length,
611
- "top_hit": {
612
- "accession": top_hit["accession"],
613
- "description": top_hit["description"],
614
- "evalue": top_hit["evalue"],
615
- "evalue_raw": str(top_hit["evalue"]),
616
- "identity_pct": top_hit["identity_pct"],
617
- "bit_score": top_hit["bit_score"],
618
- "alignment_length": top_hit.get("alignment_length", 0),
619
- } if top_hit else None,
620
- "hits": [
621
- {
622
- "accession": h["accession"],
623
- "description": h["description"],
624
- "organism": h.get("organism", ""),
625
- "evalue": h["evalue"],
626
- "evalue_raw": str(h["evalue"]),
627
- "identity_pct": h["identity_pct"],
628
- "bit_score": h["bit_score"],
629
- "alignment_length": h.get("alignment_length", 0),
630
- "query_coverage_pct": round(h.get("alignment_length", 0) / query_length * 100, 1) if query_length > 0 else 0,
631
- "hit_alignment": h.get("hit_alignment", ""),
632
- "query_alignment": h.get("query_alignment", ""),
633
- "midline": h.get("midline", ""),
634
- "score": h.get("score", 0),
635
- "positive": h.get("positive", 0),
636
- "gaps": h.get("gaps", 0),
637
- "query_from": h.get("query_from", 0),
638
- "query_to": h.get("query_to", 0),
639
- "hit_from": h.get("hit_from", 0),
640
- "hit_to": h.get("hit_to", 0),
641
- }
642
- for h in hits[:display_limit]
643
- ],
644
- }
645
-
646
-
647
- async def _run_ebi_blast_fallback(
648
- sequence: str,
649
- program: str,
650
- database: str,
651
- seq_type: str,
652
- max_hits: int,
653
- ) -> dict | None:
654
- """Run BLAST against EBI (uncached) and return a canonical result, or None
655
- if the database has no EBI equivalent or EBI itself fails."""
656
- ebi_database = EBI_BLAST_DATABASE_MAP.get(database)
657
- if not ebi_database:
658
- logger.warning("No EBI database equivalent for '%s' — skipping fallback", database)
659
- return None
660
- try:
661
- result = await BlastTool().run_uncached({
662
- "sequence": sequence,
663
- "program": program,
664
- "database": ebi_database,
665
- "max_hits": max_hits,
666
- })
667
- except Exception as e:
668
- logger.warning("EBI BLAST fallback failed: %s", e)
669
- return None
670
- if result.get("error") or not result.get("hits"):
671
- logger.warning("EBI BLAST fallback returned no hits: %s", result.get("error", "empty"))
672
- return None
673
- query_length = len("".join(sequence.replace("\n", "").replace(" ", "").split("-")))
674
- return _build_blast_result(
675
- result["hits"],
676
- source="ebi",
677
- database=database,
678
- program=program,
679
- seq_type=seq_type,
680
- query_accession="",
681
- query_length=query_length,
682
- )
683
-
684
-
685
- async def _run_blast(
686
- sequence: str,
687
- status_callback=None,
688
- fast_mode: bool = False,
689
- blast_params: dict | None = None,
690
- ) -> dict:
691
- blast_params = blast_params or {}
692
- try:
693
- program, database, seq_type = resolve_blast_params(
694
- sequence,
695
- program=blast_params.get("program"),
696
- database=blast_params.get("database"),
697
- fast_mode=fast_mode,
698
- )
699
- except ValueError as e:
700
- logger.warning("BLAST param resolution failed: %s", e)
701
- return {"error": str(e), "count": 0, "hits": []}
702
-
703
- try:
704
- max_hits = int(blast_params.get("max_hits") or 100)
705
- except (TypeError, ValueError):
706
- max_hits = 100
707
- max_hits = max(5, min(max_hits, 100))
708
- query_accession = (blast_params.get("query_accession") or "").strip()
709
-
710
- if status_callback:
711
- try:
712
- await status_callback("submitted_to_ncbi")
713
- except Exception:
714
- pass
715
-
716
- # EBI is fast & reliable (~30s) and returns UniProt accessions the rest of
717
- # the pipeline consumes directly. NCBI chronically stalls (long WAITING
718
- # even with an API key), so run EBI FIRST and fall back to NCBI.
719
- ebi_result = await _run_ebi_blast_fallback(sequence, program, database, seq_type, max_hits)
720
- if ebi_result is not None:
721
- if status_callback:
722
- try:
723
- await status_callback("parsing")
724
- except Exception:
725
- pass
726
- return ebi_result
727
-
728
- # EBI unavailable for this database or failed — try NCBI before failing.
729
- ncbi_error: str | None = None
730
- results = await ncbi_blast.run_blast_with_retry(
731
- sequence,
732
- retries=2,
733
- max_wait_seconds=600 if fast_mode else 900,
734
- database=database,
735
- program=program,
736
- hitlist_size=max_hits,
737
- )
738
-
739
- if "error" not in results:
740
- parsed = parse_blast_xml(results["raw"])
741
- if "error" not in parsed:
742
- if status_callback:
743
- try:
744
- await status_callback("parsing")
745
- except Exception:
746
- pass
747
- hits = parsed.get("hits", [])[:max_hits]
748
- return _build_blast_result(
749
- hits,
750
- source="ncbi",
751
- database=database,
752
- program=program,
753
- seq_type=seq_type,
754
- query_accession=query_accession,
755
- query_length=parsed.get("query_length", 0),
756
- )
757
- ncbi_error = parsed["error"]
758
- else:
759
- ncbi_error = results["error"]
760
-
761
- logger.warning("EBI BLAST unavailable and NCBI failed (%s)", ncbi_error)
762
- return {"error": ncbi_error or "BLAST failed via EBI and NCBI", "count": 0, "hits": []}
763
-
764
-
765
- async def _run_uniprot(top_hit: dict, query_sequence: str | None = None, try_sequence: bool = True) -> dict:
766
- accession = (top_hit.get("accession") or "").strip()
767
- if not accession:
768
- return {"error": "No accession"}
769
- desc = (top_hit.get("description") or "").strip()
770
- organism = (top_hit.get("organism") or "").strip()
771
-
772
- # Resolve ANY BLAST-hit accession (RefSeq, GenBank, PDB, Ensembl, ...) to a
773
- # UniProt accession via the strategy ladder in identifier_resolution:
774
- # direct -> xref search (fast) -> name search -> EBI sequence BLAST -> idmapping
775
- # Exhausting all tiers returns an explicit unresolved/de_novo result (§1).
776
- try:
777
- resolved = await resolve_to_uniprot(
778
- accession=accession,
779
- sequence=query_sequence,
780
- description=desc,
781
- organism=organism or None,
782
- try_sequence=try_sequence,
783
- )
784
- except Exception as e:
785
- logger.warning("Identifier resolution failed for %s: %s", accession, e)
786
- resolved = None
787
-
788
- if not resolved or resolved.get("status") != "resolved":
789
- logger.info("No UniProt mapping for %s, using BLAST data only", accession)
790
- return {
791
- "accession": accession,
792
- "full_name": top_hit.get("description", ""),
793
- "organism": top_hit.get("organism", ""),
794
- "gene_names": [],
795
- "functions": [],
796
- "keywords": [],
797
- "subcellular_locations": [],
798
- "pdb_ids": [],
799
- "go_terms": [],
800
- "sequence": "",
801
- "sequence_length": 0,
802
- "features": [],
803
- "resolution": {
804
- "uniprot_accession": None,
805
- "method": "de_novo",
806
- "original_accession": accession,
807
- },
808
- "resolved_uniprot": False,
809
- # BLAST-derived similarity is homolog-grade certainty at best (§1.1)
810
- "confidence": "homolog",
811
- "_note": f"No UniProt mapping found for {accession} — showing BLAST-derived data only",
812
- }
813
-
814
- uniprot_acc = resolved["accession"]
815
- method = resolved["method"]
816
-
817
- tool = UniprotTool()
818
- result = await tool.run({"accession": uniprot_acc})
819
- if "error" in result:
820
- return {"error": result["error"]}
821
-
822
- return {
823
- "accession": result.get("accession", uniprot_acc),
824
- "full_name": result.get("full_name", ""),
825
- "organism": result.get("organism", ""),
826
- "gene_names": result.get("gene_names", []),
827
- "functions": result.get("functions", []),
828
- "keywords": result.get("keywords", []),
829
- "subcellular_locations": result.get("subcellular_locations", []),
830
- "pdb_ids": result.get("pdb_ids", []),
831
- "go_terms": result.get("go_terms", []),
832
- "sequence": result.get("sequence", ""),
833
- "sequence_length": result.get("sequence_length", 0),
834
- "features": [
835
- f for f in (result.get("features", []) or [])
836
- if f.get("type") in ("ACTIVE_SITE", "BINDING", "MUTAGENESIS", "SITE", "MOD_RES")
837
- ],
838
- "resolution": {
839
- "uniprot_accession": uniprot_acc,
840
- "method": method,
841
- "original_accession": accession,
842
- },
843
- "resolved_uniprot": True,
844
- "confidence": resolved.get("confidence", "identified"),
845
- }
846
-
847
-
848
- async def _run_msa(query_sequence: str, blast_hits: list, alignment_mode: str = "global") -> dict:
849
- sequences = [("query", query_sequence)]
850
-
851
- for hit in blast_hits[:5]:
852
- acc = hit.get("accession", "")
853
- hit_seq = hit.get("hit_alignment", "")
854
-
855
- if acc:
856
- source = detect_source_from_accession(acc)
857
- mapped_acc = acc
858
- if source == "ncbi":
859
- mapped = await map_refseq_to_uniprot(acc)
860
- if mapped:
861
- mapped_acc = mapped
862
- try:
863
- tool = UniprotTool()
864
- ud = await tool.run({"accession": mapped_acc})
865
- if "error" not in ud and ud.get("sequence"):
866
- clean_seq = "".join(c for c in ud["sequence"] if c.isalpha()).upper()
867
- if len(clean_seq) > 10:
868
- sequences.append((acc, clean_seq))
869
- continue
870
- except Exception:
871
- pass
872
-
873
- if hit_seq:
874
- clean = "".join(c for c in hit_seq if c.isalpha()).upper()
875
- if len(clean) > 10:
876
- sequences.append((f"{acc}_aln", clean))
877
-
878
- if len(sequences) < 2:
879
- return {"error": "Not enough sequences for MSA", "aln_fasta": None, "phylotree": None}
880
-
881
- fasta_lines = []
882
- for sid, sseq in sequences:
883
- fasta_lines.append(f">{sid}")
884
- for i in range(0, len(sseq), 80):
885
- fasta_lines.append(sseq[i:i + 80])
886
- fasta_str = "\n".join(fasta_lines)
887
-
888
- try:
889
- email = settings.NCBI_EMAIL or "bioflow@example.com"
890
- seq_type = detect_sequence_type(query_sequence) or "protein"
891
- stype = "protein" if seq_type == "protein" else "dna"
892
- try:
893
- # Try local MAFFT first (fast, no network dependency)
894
- from app.tools.mafft_local import run_local_mafft
895
- local_result = await asyncio.get_running_loop().run_in_executor(
896
- None, run_local_mafft, fasta_str, "auto", 1, 300,
897
- )
898
- if local_result and local_result.get("aln_fasta"):
899
- method = "mafft-local"
900
- aln_fasta = local_result["aln_fasta"]
901
- phylotree = ""
902
- else:
903
- raise ValueError("local MAFFT unavailable or returned empty")
904
- except Exception:
905
- try:
906
- result = await run_ebi_msa_best_effort(
907
- sequence=fasta_str,
908
- stype=stype,
909
- email=email,
910
- )
911
- method = result["method"]
912
- aln_fasta = result["aln_fasta"]
913
- phylotree = result["phylotree"]
914
- except Exception as e:
915
- logger.warning("EBI MSA unavailable (%s) — using in-process fallback", e)
916
- from app.tools.msa_fallback import progressive_msa
917
- fallback = await asyncio.get_running_loop().run_in_executor(
918
- None, progressive_msa, sequences, stype
919
- )
920
- aln_fasta, phylotree = fallback
921
- method = "in-process fallback"
922
-
923
- payload = {
924
- "aln_fasta": aln_fasta,
925
- "phylotree": phylotree,
926
- "sequence_count": len(sequences),
927
- "alignment_mode": alignment_mode,
928
- "method": method,
929
- "_fallback": method != "clustalo",
930
- }
931
-
932
- # Local mode: refine query vs the best non-query sequence with an
933
- # in-process Smith-Waterman alignment so the wizard can show which
934
- # region actually matches (EBI MSA itself is always global).
935
- if alignment_mode == "local" and len(sequences) > 1:
936
- subject_id, subject_seq = sequences[1]
937
- try:
938
- payload["pairwise"] = pairwise_align(query_sequence, subject_seq, mode="local")
939
- payload["pairwise_subject"] = subject_id
940
- except Exception as e:
941
- logger.warning("Local pairwise refinement failed: %s", e)
942
-
943
- return payload
944
-
945
- except Exception as e:
946
- return {"error": str(e), "aln_fasta": None, "phylotree": None}
947
-
948
-
949
- async def _run_domains(accession: str) -> dict:
950
- """Run domain analysis using the shared tool module (eliminates code duplication)."""
951
- try:
952
- from app.tools.domain_analysis import fetch_interpro_domains
953
- return await fetch_interpro_domains(accession)
954
- except Exception as e:
955
- return {"error": str(e), "uniprot_accession": accession, "sequence_length": 0, "domains": []}
956
-
957
-
958
- async def _run_denovo_steps(job_id: str, sequence: str, steps: list[str], _mark, context: dict) -> None:
959
- """Tier-6 branch (techspec.md §1): characterize from sequence alone.
960
-
961
- Runs when BLAST finds no homolog at all. Composition/function hints land
962
- in the uniprot slot (the annotation slot); MSA/phylo/pathways are marked
963
- explicitly unavailable rather than left empty.
964
- """
965
- import asyncio as _asyncio
966
- from app.services.de_novo import (
967
- composition_stats, esmfold_structure, function_hints, interpro_sequence_search,
968
- )
969
-
970
- unavailable = "Unavailable for de novo sequences — no identified homolog"
971
-
972
- async def _do_annotation():
973
- bundle = {"_de_novo": True}
974
- try:
975
- bundle["composition"] = composition_stats(sequence)
976
- except Exception as e:
977
- bundle["composition"] = {"error": str(e)}
978
- try:
979
- bundle["function_hints"] = function_hints(sequence)
980
- except Exception as e:
981
- bundle["function_hints"] = {"error": str(e)}
982
- return bundle
983
-
984
- async def _do_domains():
985
- try:
986
- return await interpro_sequence_search(sequence)
987
- except Exception as e:
988
- logger.warning("De novo InterProScan failed: %s", e)
989
- return {"error": str(e), "domains": [], "sequence_length": 0}
990
-
991
- async def _do_structure():
992
- try:
993
- return await esmfold_structure(sequence)
994
- except Exception as e:
995
- logger.warning("De novo ESMFold failed: %s", e)
996
- return {"structure_available": False, "source": "esmfold", "message": str(e)}
997
-
998
- fan = []
999
- names = []
1000
- if "uniprot" in steps:
1001
- fan.append(_do_annotation())
1002
- names.append("uniprot")
1003
- if "domains" in steps:
1004
- fan.append(_do_domains())
1005
- names.append("domains")
1006
- if "alphafold" in steps:
1007
- fan.append(_do_structure())
1008
- names.append("alphafold")
1009
-
1010
- for name in names:
1011
- _mark(name, "running", progress=10)
1012
- results = await _asyncio.gather(*fan, return_exceptions=True)
1013
-
1014
- for name, res in zip(names, results):
1015
- if isinstance(res, Exception):
1016
- _mark(name, "failed", error=str(res)[:500])
1017
- continue
1018
- ok = (
1019
- res.get("domains") is not None if name == "domains"
1020
- else res.get("structure_available") is True if name == "alphafold"
1021
- else True
1022
- )
1023
- _mark(name, "complete" if ok else "failed", progress=100,
1024
- data=res, error=None if ok else "Prediction returned nothing usable")
1025
- context[name] = res
1026
-
1027
- # Annotation-database features have no de novo substitute — say so.
1028
- for name in ("msa", "phylo", "pathway_enrichment"):
1029
- if name in steps:
1030
- _mark(name, "failed", error=unavailable)
1031
-
1032
-
1033
- async def _run_domains_or_denovo(sequence: str, accession: str | None, resolved_uniprot: bool) -> dict:
1034
- """Accession lookup when resolved; InterProScan sequence-search otherwise."""
1035
- if accession and resolved_uniprot:
1036
- return await _run_domains(accession)
1037
- from app.services.de_novo import interpro_sequence_search
1038
- try:
1039
- result = await interpro_sequence_search(sequence)
1040
- result["confidence"] = "de_novo"
1041
- return result
1042
- except Exception as e:
1043
- return {"error": str(e), "domains": [], "sequence_length": 0}
1044
-
1045
-
1046
- async def _run_alphafold_or_esmfold(
1047
- context: dict, sequence: str, accession: str | None, resolved_uniprot: bool,
1048
- ) -> dict:
1049
- """AlphaFold DB lookup when resolved; ESMFold ab initio otherwise."""
1050
- if accession and resolved_uniprot:
1051
- result = await _run_alphafold(context)
1052
- return result or {}
1053
- from app.services.de_novo import esmfold_structure
1054
- try:
1055
- return await esmfold_structure(sequence)
1056
- except Exception as e:
1057
- return {"structure_available": False, "source": "esmfold", "message": str(e)}
1058
-
1059
-
1060
- async def _run_pathway_enrichment(context: dict) -> dict | None:
1061
- gene_names = []
1062
- uniprot = context.get("uniprot", {})
1063
- if isinstance(uniprot, dict):
1064
- gene_names = uniprot.get("gene_names", [])[:20] if isinstance(uniprot.get("gene_names"), list) else []
1065
- if not gene_names:
1066
- blast_data = context.get("blast", {})
1067
- if isinstance(blast_data, dict):
1068
- for hit in (blast_data.get("hits") or [])[:10]:
1069
- words = (hit.get("description", "") or "").replace("(", " ").replace(")", " ").split()
1070
- for w in words:
1071
- if w.isupper() and len(w) >= 2 and not w.startswith("OS="):
1072
- gene_names.append(w)
1073
- break
1074
- if not gene_names:
1075
- return None
1076
- try:
1077
- from app.services.pathway_enrichment import run_enrichment
1078
- result = await run_enrichment(gene_names)
1079
- return result
1080
- except Exception as e:
1081
- logger.warning(f"Pathway enrichment failed: {e}")
1082
- return None
1083
-
1084
-
1085
- async def _run_alphafold(context: dict) -> dict | None:
1086
- uniprot_data = context.get("uniprot", {})
1087
- accession = uniprot_data.get("accession") if isinstance(uniprot_data, dict) else None
1088
- if not accession:
1089
- return None
1090
- try:
1091
- from app.tools.alphafold import AlphaFoldTool
1092
- result = await AlphaFoldTool().run({"uniprot_accession": accession})
1093
- return result
1094
- except Exception as e:
1095
- logger.warning(f"AlphaFold fetch failed for {accession}: {e}")
1096
- return {"structure_available": False, "message": str(e)}
1097
-
1098
-
1099
- async def _run_interpret(context: dict) -> dict:
1100
- from app.ai.interpreter import interpret_text
1101
- prompt_context = {
1102
- "blast": context.get("blast", {}),
1103
- "uniprot": context.get("uniprot", {}),
1104
- "alphafold": context.get("alphafold", {}),
1105
- "pathway_enrichment": context.get("pathway_enrichment", {}),
1106
- "query_confidence": context.get("query", {}).get("confidence", "identified"),
1107
- }
1108
- return await interpret_text("protein_analysis", prompt_context)
 
1
+ """
2
+ In-memory pipeline v2 — runs BLAST → UniProt → MSA → Phylo → Domains → Interpretation
3
+ in a background thread. Uses a thread-safe dict for job storage.
4
+ """
5
+
6
+ import asyncio
7
+ import logging
8
+ import threading
9
+ import uuid
10
+ from datetime import datetime, timezone
11
+
12
+ from fastapi import APIRouter, HTTPException, Depends, Request
13
+ from pydantic import BaseModel, Field
14
+
15
+ from app.config import settings
16
+ from app.deps import limiter
17
+ from app.services.rate_limit import check_daily_limit_pipelines
18
+ from app.services.auth import get_user_id
19
+ from app.services.supabase import get_supabase
20
+ from app.integrations.ncbi import blast as ncbi_blast
21
+ from app.integrations.ncbi.parser import parse_blast_xml
22
+ from app.services.validators import validate_fasta
23
+ from app.services.sequence_utils import detect_source_from_accession, map_refseq_to_uniprot, detect_sequence_type
24
+ from app.services.blast_config import resolve_blast_params
25
+ from app.services.identifier_resolution import resolve_to_uniprot
26
+ from app.tools.ebi_msa import EBI_TOOLS, run_ebi_msa, run_ebi_msa_best_effort
27
+ from app.tools.blast import BlastTool
28
+ from app.tools.uniprot import UniprotTool
29
+ from app.tools.pairwise_alignment import pairwise_align, VALID_MODES
30
+
31
+ logger = logging.getLogger(__name__)
32
+ router = APIRouter()
33
+
34
+ _jobs: dict[str, dict] = {}
35
+ _jobs_lock = threading.Lock()
36
+
37
+ STEP_ORDER = ["blast", "uniprot", "msa", "phylo", "domains", "pathway_enrichment", "alphafold", "interpret"]
38
+
39
+
40
+ def _get_job(job_id: str) -> dict | None:
41
+ with _jobs_lock:
42
+ return _jobs.get(job_id)
43
+
44
+
45
+ def _is_real_uuid(job_id: str) -> bool:
46
+ try:
47
+ uuid.UUID(job_id)
48
+ return True
49
+ except (ValueError, AttributeError):
50
+ return False
51
+
52
+
53
+ def _persist_v2_final(job_id: str, status: str, context: dict, error: str | None = None):
54
+ """Update the persisted wizard job row with the final status + full context."""
55
+ if not _is_real_uuid(job_id):
56
+ return
57
+ now = datetime.now(timezone.utc).isoformat()
58
+ payload = {
59
+ "status": status,
60
+ "progress_pct": 100,
61
+ "completed_at": now,
62
+ "context_json": context,
63
+ "error_message": error,
64
+ }
65
+ if error:
66
+ payload["error"] = error
67
+ try:
68
+ get_supabase().table("jobs").update(payload).eq("id", job_id).execute()
69
+ except Exception as e:
70
+ logger.warning("Could not update v2 job %s to Supabase: %s", job_id, e)
71
+
72
+
73
+ def _persist_v2_job(job_id: str, payload: dict):
74
+ """Persist a wizard v2 job to Supabase (best-effort) so it shows up in
75
+ the user's job history and can be shared. Uses the service role key."""
76
+ try:
77
+ get_supabase().table("jobs").insert(payload).execute()
78
+ except Exception as e:
79
+ logger.warning("Could not persist v2 job %s to Supabase: %s", job_id, e)
80
+
81
+
82
+ def _set_step_status(job_id: str, step: str, status: str, progress: int = 0, data: dict | None = None, error: str | None = None):
83
+ with _jobs_lock:
84
+ if job_id not in _jobs:
85
+ return
86
+ _jobs[job_id]["steps"][step] = {"status": status, "progress": progress, "data": data, "error": error}
87
+ if status == "running":
88
+ _jobs[job_id]["current_step"] = step
89
+
90
+
91
+ def _set_job_failed(job_id: str, message: str):
92
+ with _jobs_lock:
93
+ if job_id in _jobs:
94
+ _jobs[job_id]["status"] = "failed"
95
+ _jobs[job_id]["error"] = message
96
+
97
+
98
+ class PipelineV2RunRequest(BaseModel):
99
+ sequence: str = Field(..., min_length=6, description="Protein sequence (FASTA or raw)")
100
+ steps: list[str] = Field(default_factory=lambda: list(STEP_ORDER), description="Steps to run")
101
+ fast_mode: bool = Field(default=False, description="Use Swiss-Prot instead of nr for faster results")
102
+ database: str = Field("", description="BLAST database override")
103
+ program: str = Field("", description="BLAST program override")
104
+ max_hits: int = Field(100, description="Max BLAST hits to return")
105
+ query_accession: str = Field("", description="Optional query accession for display")
106
+ alignment_mode: str = Field("global", description="Alignment mode for the MSA step: 'global' (full-length) or 'local' (Smith-Waterman refinement of query vs top hit)")
107
+ parent_job_id: str | None = Field(None, description="Parent job ID for DAG branching")
108
+
109
+
110
+ @router.post("/run")
111
+ async def run_pipeline_v2(request: Request, req: PipelineV2RunRequest, user_id: str | None = Depends(get_user_id)):
112
+ validation = validate_fasta(req.sequence, "blast")
113
+ if not validation.valid:
114
+ raise HTTPException(status_code=400, detail=validation.error)
115
+
116
+ alignment_mode = req.alignment_mode.lower()
117
+ if alignment_mode not in VALID_MODES:
118
+ raise HTTPException(status_code=400, detail=f"alignment_mode must be one of {list(VALID_MODES)}, got {alignment_mode!r}")
119
+
120
+ seq = str(validation.sequences[0].seq).upper()
121
+ clean = "".join(c for c in seq if c.isalpha())
122
+
123
+ job_id = str(uuid.uuid4())
124
+ now = datetime.now(timezone.utc).isoformat()
125
+
126
+ requested = [s for s in req.steps if s in STEP_ORDER]
127
+ if not requested:
128
+ requested = list(STEP_ORDER)
129
+
130
+ steps_dict = {s: {"status": "pending", "progress": 0, "data": None, "error": None} for s in STEP_ORDER}
131
+
132
+ blast_params = {
133
+ "database": req.database,
134
+ "program": req.program,
135
+ "max_hits": req.max_hits,
136
+ "query_accession": req.query_accession,
137
+ }
138
+
139
+ with _jobs_lock:
140
+ _jobs[job_id] = {
141
+ "job_id": job_id,
142
+ "status": "running",
143
+ "current_step": None,
144
+ "steps": steps_dict,
145
+ "requested_steps": requested,
146
+ "sequence": clean,
147
+ "blast_params": blast_params,
148
+ "alignment_mode": alignment_mode,
149
+ "error": None,
150
+ "created_at": now,
151
+ }
152
+
153
+ # Best-effort persistence so wizard jobs appear in history and can be shared.
154
+ persist_payload = {
155
+ "id": job_id,
156
+ "user_id": user_id,
157
+ "tool": "wizard_v2",
158
+ "pipeline_type": "wizard_v2",
159
+ "query_preview": clean[:60],
160
+ "status": "running",
161
+ "progress_pct": 0,
162
+ "title": "Wizard pipeline",
163
+ "description": "Pipeline v2 wizard run",
164
+ "created_at": now,
165
+ }
166
+ if req.parent_job_id:
167
+ persist_payload["parent_job_id"] = req.parent_job_id
168
+ _persist_v2_job(job_id, persist_payload)
169
+
170
+ t = threading.Thread(
171
+ target=_run_pipeline,
172
+ args=(job_id, clean, requested),
173
+ kwargs={"fast_mode": req.fast_mode, "blast_params": blast_params, "alignment_mode": alignment_mode},
174
+ daemon=True,
175
+ )
176
+ t.start()
177
+
178
+ return {"job_id": job_id}
179
+
180
+
181
+ @router.get("/status/{job_id}")
182
+ async def get_pipeline_v2_status(job_id: str):
183
+ job = _get_job(job_id)
184
+ if not job:
185
+ raise HTTPException(status_code=404, detail="Job not found")
186
+ return job
187
+
188
+
189
+ async def run_pipeline(
190
+ sequence: str,
191
+ organism: str = "Homo sapiens",
192
+ analysis_type: str = "comprehensive",
193
+ status_callback=None,
194
+ fast_mode: bool = False,
195
+ blast_params: dict | None = None,
196
+ ) -> dict:
197
+ """Public async entry point for the pipeline (used by pipeline_worker).
198
+
199
+ Creates a temporary in-memory job, runs the configured steps, and
200
+ returns the context dict with all results.
201
+ """
202
+ job_id = f"worker-{uuid.uuid4().hex[:12]}"
203
+ requested = list(STEP_ORDER)
204
+ steps_dict = {s: {"status": "pending", "progress": 0, "data": None, "error": None} for s in STEP_ORDER}
205
+
206
+ with _jobs_lock:
207
+ _jobs[job_id] = {
208
+ "job_id": job_id,
209
+ "status": "running",
210
+ "current_step": None,
211
+ "steps": steps_dict,
212
+ "requested_steps": requested,
213
+ "sequence": sequence,
214
+ "blast_params": blast_params or {},
215
+ "error": None,
216
+ "created_at": datetime.now(timezone.utc).isoformat(),
217
+ }
218
+
219
+ try:
220
+ await _execute(
221
+ job_id,
222
+ sequence,
223
+ requested,
224
+ status_callback=status_callback,
225
+ fast_mode=fast_mode,
226
+ blast_params=blast_params,
227
+ )
228
+ finally:
229
+ job = _get_job(job_id)
230
+ with _jobs_lock:
231
+ _jobs.pop(job_id, None)
232
+
233
+ if job and job.get("status") == "failed":
234
+ raise RuntimeError(job.get("error", "Pipeline failed"))
235
+
236
+ query_accession = ((blast_params or {}).get("query_accession") or "").strip()
237
+ context: dict = {
238
+ "sequence": sequence,
239
+ "length": len(sequence),
240
+ "query": {
241
+ "sequence": sequence,
242
+ "length": len(sequence),
243
+ "sequence_type": detect_sequence_type(sequence) or "protein",
244
+ },
245
+ }
246
+ if query_accession:
247
+ context["query"]["accession"] = query_accession
248
+ if job:
249
+ for step_name, step_info in job.get("steps", {}).items():
250
+ if step_info.get("data"):
251
+ context[step_name] = step_info["data"]
252
+ return context
253
+
254
+
255
+ # ---------------------------------------------------------------------------
256
+ # Background pipeline
257
+ # ---------------------------------------------------------------------------
258
+
259
+ def _run_pipeline(job_id: str, sequence: str, steps: list[str], fast_mode: bool = False, blast_params: dict | None = None, alignment_mode: str = "global"):
260
+ loop = asyncio.new_event_loop()
261
+ asyncio.set_event_loop(loop)
262
+ try:
263
+ loop.run_until_complete(
264
+ _execute(job_id, sequence, steps, fast_mode=fast_mode, blast_params=blast_params, alignment_mode=alignment_mode)
265
+ )
266
+ except Exception as e:
267
+ logger.exception(f"[{job_id}] Unhandled pipeline error")
268
+ _set_job_failed(job_id, f"Pipeline error: {e}")
269
+ finally:
270
+ loop.close()
271
+ asyncio.set_event_loop(None)
272
+
273
+
274
+ async def _execute(job_id: str, sequence: str, steps: list[str], status_callback=None, fast_mode: bool = False, blast_params: dict | None = None, alignment_mode: str = "global"):
275
+ context: dict = {
276
+ "sequence": sequence,
277
+ "length": len(sequence),
278
+ "query": {
279
+ "sequence": sequence,
280
+ "length": len(sequence),
281
+ "sequence_type": detect_sequence_type(sequence) or "protein",
282
+ },
283
+ }
284
+
285
+ _STEP_FRONTEND = {
286
+ "blast": "running",
287
+ "uniprot": "fetching_uniprot",
288
+ "msa": "running_msa",
289
+ "phylo": "running_msa",
290
+ "domains": "fetching_uniprot",
291
+ "pathway_enrichment": "pathway_enrichment",
292
+ "alphafold": "fetching_alphafold",
293
+ "interpret": "interpreting",
294
+ }
295
+
296
+ _failed_step = None
297
+ _failed_error = None
298
+
299
+ async def _notify(step_key: str):
300
+ if status_callback:
301
+ try:
302
+ await status_callback(_STEP_FRONTEND.get(step_key, "running"))
303
+ except Exception:
304
+ pass
305
+
306
+ def _mark(step_key: str, status: str, **kw):
307
+ _set_step_status(job_id, step_key, status, **kw)
308
+
309
+ def _fail(step_key: str, msg: str):
310
+ nonlocal _failed_step, _failed_error
311
+ _mark(step_key, "failed", error=msg)
312
+ _failed_step = step_key
313
+ _failed_error = msg
314
+
315
+ # ---- Step 1: BLAST (must run first) ----
316
+ denovo_mode = False
317
+ if "blast" in steps:
318
+ await _notify("blast")
319
+ _mark("blast", "running", progress=10)
320
+ result = await _run_blast(
321
+ sequence,
322
+ status_callback=status_callback,
323
+ fast_mode=fast_mode,
324
+ blast_params=blast_params,
325
+ )
326
+ zero_hits = result.get("count", 0) == 0
327
+ if zero_hits and (detect_sequence_type(sequence) or "protein") == "protein":
328
+ # Tier 6: no database match at all — characterize from sequence
329
+ # alone instead of failing the run (techspec.md §1).
330
+ denovo_mode = True
331
+ result["_note"] = "No BLAST hits found — switching to de novo characterization"
332
+ _mark("blast", "complete", progress=100, data=result)
333
+ else:
334
+ _mark("blast", "failed" if zero_hits else "complete", progress=100, data=result)
335
+ if zero_hits:
336
+ _failed_step = "blast"
337
+ _failed_error = result.get("error", "No BLAST hits found")
338
+ context["blast"] = result
339
+
340
+ if denovo_mode:
341
+ context["query"]["confidence"] = "de_novo"
342
+ await _run_denovo_steps(job_id, sequence, steps, _mark, context)
343
+ await _finalize_context(job_id, context)
344
+ with _jobs_lock:
345
+ if job_id in _jobs:
346
+ _jobs[job_id]["status"] = "complete"
347
+ _jobs[job_id]["context"] = context
348
+ _persist_v2_final(job_id, "complete", context)
349
+ return
350
+
351
+ # ---- Step 2: Fan-out — UniProt, MSA, Pathway run in parallel ----
352
+ # They all only depend on BLAST results, not on each other.
353
+ blast_data = context.get("blast", {})
354
+ hits = (blast_data.get("hits") if isinstance(blast_data, dict) else []) or []
355
+ top_hit = blast_data.get("top_hit") if isinstance(blast_data, dict) else None
356
+
357
+ async def _do_uniprot():
358
+ candidates = ([top_hit] + hits[:5]) if top_hit else hits[:5]
359
+ last_result: dict | None = None
360
+ for idx, candidate in enumerate(candidates):
361
+ # Only the top candidate may trigger the expensive EBI sequence
362
+ # BLAST fallback — otherwise a pathological query could run it up
363
+ # to 6 times (up to ~3 min each).
364
+ last_result = await _run_uniprot(
365
+ candidate,
366
+ query_sequence=sequence,
367
+ try_sequence=idx == 0,
368
+ )
369
+ if "error" not in last_result:
370
+ return last_result
371
+ return last_result if last_result else {"error": "No BLAST hits for UniProt lookup"}
372
+
373
+ async def _do_msa():
374
+ if not hits:
375
+ return {"error": "No BLAST hits for MSA"}
376
+ return await _run_msa(sequence, hits, alignment_mode)
377
+
378
+ async def _do_pathway():
379
+ return await _run_pathway_enrichment(context)
380
+
381
+ fan_out = []
382
+ fan_names = []
383
+ if "uniprot" in steps and not _failed_step:
384
+ fan_out.append(_do_uniprot())
385
+ fan_names.append("uniprot")
386
+ if "msa" in steps and not _failed_step:
387
+ fan_out.append(_do_msa())
388
+ fan_names.append("msa")
389
+ if "pathway_enrichment" in steps and not _failed_step:
390
+ fan_out.append(_do_pathway())
391
+ fan_names.append("pathway_enrichment")
392
+
393
+ if fan_out:
394
+ # Notify for the first active step in the fan-out
395
+ await _notify(fan_names[0])
396
+ for name in fan_names:
397
+ _mark(name, "running", progress=10)
398
+
399
+ results = await asyncio.gather(*fan_out, return_exceptions=True)
400
+
401
+ for name, res in zip(fan_names, results):
402
+ if isinstance(res, Exception):
403
+ _fail(name, str(res)[:500])
404
+ continue
405
+
406
+ if name == "uniprot":
407
+ s = "complete" if "error" not in res else "failed"
408
+ _mark("uniprot", s, progress=100, data=res)
409
+ context["uniprot"] = res
410
+ if res.get("confidence"):
411
+ context["query"]["confidence"] = res["confidence"]
412
+ if "error" in res:
413
+ _failed_step = "uniprot"
414
+ _failed_error = res["error"]
415
+
416
+ elif name == "msa":
417
+ s = "complete" if res.get("aln_fasta") else "failed"
418
+ _mark("msa", s, progress=100, data=res)
419
+ context["msa"] = res
420
+ if res.get("phylotree"):
421
+ context.setdefault("phylo_data", {})["phylotree_newick"] = res["phylotree"]
422
+
423
+ elif name == "pathway_enrichment":
424
+ s = "complete" if res and res.get("pathways") else "failed"
425
+ _mark("pathway_enrichment", s, progress=100, data=res or {})
426
+ context["pathway_enrichment"] = res
427
+
428
+ # ---- Step 3: Phylo (instant — copies from MSA) ----
429
+ if "phylo" in steps and not _failed_step:
430
+ _mark("phylo", "running", progress=10)
431
+ newick = None
432
+ msa_data = context.get("msa", {})
433
+ if isinstance(msa_data, dict):
434
+ newick = msa_data.get("phylotree")
435
+ if not newick:
436
+ newick = context.get("phylo_data", {}).get("phylotree_newick")
437
+ if newick:
438
+ _mark("phylo", "complete", progress=100, data={"phylotree_newick": newick})
439
+ context["phylo"] = {"phylotree_newick": newick}
440
+ else:
441
+ _mark("phylo", "failed", error="No phylotree available from MSA")
442
+
443
+ # ---- Step 4: Domains + AlphaFold in parallel ----
444
+ # When UniProt resolution exhausted all tiers, these fall back to their
445
+ # sequence-only equivalents instead of being silently skipped (§1.2).
446
+ uniprot_data = context.get("uniprot", {})
447
+ accession = uniprot_data.get("accession") if isinstance(uniprot_data, dict) else None
448
+ resolved_uniprot = bool(uniprot_data.get("resolved_uniprot")) if isinstance(uniprot_data, dict) else False
449
+
450
+ post_uniprot = []
451
+ post_uniprot_names = []
452
+ if "domains" in steps and not _failed_step:
453
+ post_uniprot.append(
454
+ _run_domains_or_denovo(sequence, accession, resolved_uniprot)
455
+ )
456
+ post_uniprot_names.append("domains")
457
+ if "alphafold" in steps and not _failed_step:
458
+ post_uniprot.append(
459
+ _run_alphafold_or_esmfold(context, sequence, accession, resolved_uniprot)
460
+ )
461
+ post_uniprot_names.append("alphafold")
462
+
463
+ if post_uniprot:
464
+ await _notify(post_uniprot_names[0])
465
+ for name in post_uniprot_names:
466
+ _mark(name, "running", progress=10)
467
+
468
+ results2 = await asyncio.gather(*post_uniprot, return_exceptions=True)
469
+
470
+ for name, res in zip(post_uniprot_names, results2):
471
+ if isinstance(res, Exception):
472
+ _fail(name, str(res)[:500])
473
+ continue
474
+
475
+ if name == "domains":
476
+ s = "complete" if res.get("domains") is not None else "failed"
477
+ _mark("domains", s, progress=100, data=res)
478
+ context["domains"] = res
479
+ elif name == "alphafold":
480
+ s = "complete" if res else "failed"
481
+ _mark("alphafold", s, progress=100, data=res or {})
482
+ context["alphafold"] = res
483
+
484
+ # ---- Step 5: Interpret (needs all context) ----
485
+ if "interpret" in steps and not _failed_step:
486
+ await _notify("interpret")
487
+ _mark("interpret", "running", progress=10)
488
+ result = await _run_interpret(context)
489
+ s = "complete" if result.get("interpretation") else "failed"
490
+ _mark("interpret", s, progress=100, data=result)
491
+ context["interpret"] = result
492
+
493
+ # ---- Final status ----
494
+ if _failed_step and _failed_step in ("blast", "uniprot"):
495
+ with _jobs_lock:
496
+ if job_id in _jobs:
497
+ _jobs[job_id]["status"] = "failed"
498
+ _jobs[job_id]["error"] = f"Pipeline failed at {_failed_step}: {_failed_error}"
499
+ _persist_v2_final(job_id, "failed", context, error=f"Pipeline failed at {_failed_step}: {_failed_error}")
500
+ else:
501
+ await _finalize_context(job_id, context)
502
+ with _jobs_lock:
503
+ if job_id in _jobs:
504
+ _jobs[job_id]["status"] = "complete"
505
+ _jobs[job_id]["context"] = context
506
+ _persist_v2_final(job_id, "complete", context)
507
+
508
+
509
+ async def _finalize_context(job_id: str, context: dict):
510
+ """Pre-persist hook: capture source pages (§3) + build the final report."""
511
+ user_id = None
512
+ with _jobs_lock:
513
+ entry = _jobs.get(job_id)
514
+ if entry:
515
+ user_id = entry.get("user_id")
516
+ try:
517
+ _capture_run_sources(job_id, context, user_id)
518
+ except Exception as e: # captures must never break a run
519
+ logger.warning("[%s] page capture wiring failed: %s", job_id, e)
520
+ try:
521
+ from app.services.final_synthesis import synthesize
522
+
523
+ context["final_report"] = await synthesize(context)
524
+ except Exception as e:
525
+ logger.warning("[%s] final synthesis failed: %s", job_id, e)
526
+
527
+
528
+ def _capture_run_sources(job_id: str, context: dict, user_id: str | None):
529
+ """Queue one page_captures row per external source this run queried."""
530
+ from app.services.page_capture import capture_bg
531
+
532
+ blast = context.get("blast") or {}
533
+ top_hit = blast.get("top_hit")
534
+ if top_hit and top_hit.get("accession"):
535
+ capture_bg(job_id, "ncbi", f"https://www.ncbi.nlm.nih.gov/protein/{top_hit['accession']}", user_id)
536
+
537
+ uniprot = context.get("uniprot") or {}
538
+ if uniprot.get("_de_novo"):
539
+ pass # no external annotation source was queried in de novo mode
540
+ elif uniprot.get("accession"):
541
+ capture_bg(job_id, "uniprot", f"https://www.uniprot.org/uniprotkb/{uniprot['accession']}", user_id)
542
+ pdb_ids = uniprot.get("pdb_ids") or []
543
+ if pdb_ids:
544
+ capture_bg(job_id, "rcsb", f"https://www.rcsb.org/structure/{pdb_ids[0]}", user_id)
545
+
546
+ domains = context.get("domains") or {}
547
+ dom_list = domains.get("domains") or []
548
+ if dom_list and dom_list[0].get("accession"):
549
+ acc = dom_list[0]["accession"]
550
+ if str(acc).startswith("IPR"):
551
+ capture_bg(job_id, "interpro", f"https://www.ebi.ac.uk/interpro/entry/InterPro/{acc}", user_id)
552
+
553
+ af = context.get("alphafold") or {}
554
+ if af.get("structure_available") and af.get("source") != "esmfold" and uniprot.get("accession"):
555
+ capture_bg(job_id, "alphafold", f"https://alphafold.ebi.ac.uk/uniprot/{uniprot['accession']}", user_id)
556
+
557
+ pathway = context.get("pathway_enrichment") or {}
558
+ pw_list = (pathway.get("pathways") if isinstance(pathway, dict) else None) or []
559
+ if pw_list and pw_list[0].get("stId"):
560
+ # run_enrichment queries the Reactome projection API exclusively
561
+ capture_bg(
562
+ job_id, "reactome",
563
+ f"https://reactome.org/PathwayBrowser/#{pw_list[0]['stId']}",
564
+ user_id,
565
+ )
566
+
567
+
568
+ # ---------------------------------------------------------------------------
569
+ # Step implementations
570
+ # ---------------------------------------------------------------------------
571
+
572
+ # NCBI database name → EBI ncbiblast database name. EBI has no direct "nr"
573
+ # equivalent; uniprotkb is its closest comprehensive protein database.
574
+ EBI_BLAST_DATABASE_MAP = {
575
+ "nr": "uniprotkb",
576
+ "swissprot": "uniprotkb_swissprot",
577
+ "pdb": "pdb",
578
+ "pdbaa": "pdb_nr",
579
+ "refseq_protein": "refseq_protein",
580
+ "env_nr": "env_nr",
581
+ "nt": "nt",
582
+ "refseq_rna": "refseq_rna",
583
+ "refseq_genomic": "refseq_genomic",
584
+ "est": "est",
585
+ "gss": "gss",
586
+ }
587
+
588
+
589
+ def _build_blast_result(
590
+ hits: list[dict],
591
+ *,
592
+ source: str,
593
+ database: str,
594
+ program: str,
595
+ seq_type: str,
596
+ query_accession: str,
597
+ query_length: int,
598
+ display_limit: int = 20,
599
+ ) -> dict:
600
+ """Normalize parsed hits (NCBI XML parser or EBI tool) into the pipeline's
601
+ canonical BLAST result shape consumed by the frontend."""
602
+ top_hit = hits[0] if hits else None
603
+ return {
604
+ "count": len(hits),
605
+ "source": source,
606
+ "database": database,
607
+ "program": program,
608
+ "query_sequence_type": seq_type,
609
+ "query_accession": query_accession,
610
+ "query_length": query_length,
611
+ "top_hit": {
612
+ "accession": top_hit["accession"],
613
+ "description": top_hit["description"],
614
+ "evalue": top_hit["evalue"],
615
+ "evalue_raw": str(top_hit["evalue"]),
616
+ "identity_pct": top_hit["identity_pct"],
617
+ "bit_score": top_hit["bit_score"],
618
+ "alignment_length": top_hit.get("alignment_length", 0),
619
+ } if top_hit else None,
620
+ "hits": [
621
+ {
622
+ "accession": h["accession"],
623
+ "description": h["description"],
624
+ "organism": h.get("organism", ""),
625
+ "evalue": h["evalue"],
626
+ "evalue_raw": str(h["evalue"]),
627
+ "identity_pct": h["identity_pct"],
628
+ "bit_score": h["bit_score"],
629
+ "alignment_length": h.get("alignment_length", 0),
630
+ "query_coverage_pct": round(h.get("alignment_length", 0) / query_length * 100, 1) if query_length > 0 else 0,
631
+ "hit_alignment": h.get("hit_alignment", ""),
632
+ "query_alignment": h.get("query_alignment", ""),
633
+ "midline": h.get("midline", ""),
634
+ "score": h.get("score", 0),
635
+ "positive": h.get("positive", 0),
636
+ "gaps": h.get("gaps", 0),
637
+ "query_from": h.get("query_from", 0),
638
+ "query_to": h.get("query_to", 0),
639
+ "hit_from": h.get("hit_from", 0),
640
+ "hit_to": h.get("hit_to", 0),
641
+ }
642
+ for h in hits[:display_limit]
643
+ ],
644
+ }
645
+
646
+
647
+ async def _run_ebi_blast_fallback(
648
+ sequence: str,
649
+ program: str,
650
+ database: str,
651
+ seq_type: str,
652
+ max_hits: int,
653
+ ) -> dict | None:
654
+ """Run BLAST against EBI (uncached) and return a canonical result, or None
655
+ if the database has no EBI equivalent or EBI itself fails."""
656
+ ebi_database = EBI_BLAST_DATABASE_MAP.get(database)
657
+ if not ebi_database:
658
+ logger.warning("No EBI database equivalent for '%s' — skipping fallback", database)
659
+ return None
660
+ try:
661
+ result = await BlastTool().run_uncached({
662
+ "sequence": sequence,
663
+ "program": program,
664
+ "database": ebi_database,
665
+ "max_hits": max_hits,
666
+ })
667
+ except Exception as e:
668
+ logger.warning("EBI BLAST fallback failed: %s", e)
669
+ return None
670
+ if result.get("error") or not result.get("hits"):
671
+ logger.warning("EBI BLAST fallback returned no hits: %s", result.get("error", "empty"))
672
+ return None
673
+ query_length = len("".join(sequence.replace("\n", "").replace(" ", "").split("-")))
674
+ return _build_blast_result(
675
+ result["hits"],
676
+ source="ebi",
677
+ database=database,
678
+ program=program,
679
+ seq_type=seq_type,
680
+ query_accession="",
681
+ query_length=query_length,
682
+ )
683
+
684
+
685
+ async def _run_blast(
686
+ sequence: str,
687
+ status_callback=None,
688
+ fast_mode: bool = False,
689
+ blast_params: dict | None = None,
690
+ ) -> dict:
691
+ blast_params = blast_params or {}
692
+ try:
693
+ program, database, seq_type = resolve_blast_params(
694
+ sequence,
695
+ program=blast_params.get("program"),
696
+ database=blast_params.get("database"),
697
+ fast_mode=fast_mode,
698
+ )
699
+ except ValueError as e:
700
+ logger.warning("BLAST param resolution failed: %s", e)
701
+ return {"error": str(e), "count": 0, "hits": []}
702
+
703
+ try:
704
+ max_hits = int(blast_params.get("max_hits") or 100)
705
+ except (TypeError, ValueError):
706
+ max_hits = 100
707
+ max_hits = max(5, min(max_hits, 100))
708
+ query_accession = (blast_params.get("query_accession") or "").strip()
709
+
710
+ if status_callback:
711
+ try:
712
+ await status_callback("submitted_to_ncbi")
713
+ except Exception:
714
+ pass
715
+
716
+ # EBI is fast & reliable (~30s) and returns UniProt accessions the rest of
717
+ # the pipeline consumes directly. NCBI chronically stalls (long WAITING
718
+ # even with an API key), so run EBI FIRST and fall back to NCBI.
719
+ ebi_result = await _run_ebi_blast_fallback(sequence, program, database, seq_type, max_hits)
720
+ if ebi_result is not None:
721
+ if status_callback:
722
+ try:
723
+ await status_callback("parsing")
724
+ except Exception:
725
+ pass
726
+ return ebi_result
727
+
728
+ # EBI unavailable for this database or failed — try NCBI before failing.
729
+ ncbi_error: str | None = None
730
+ results = await ncbi_blast.run_blast_with_retry(
731
+ sequence,
732
+ retries=2,
733
+ max_wait_seconds=600 if fast_mode else 900,
734
+ database=database,
735
+ program=program,
736
+ hitlist_size=max_hits,
737
+ )
738
+
739
+ if "error" not in results:
740
+ parsed = parse_blast_xml(results["raw"])
741
+ if "error" not in parsed:
742
+ if status_callback:
743
+ try:
744
+ await status_callback("parsing")
745
+ except Exception:
746
+ pass
747
+ hits = parsed.get("hits", [])[:max_hits]
748
+ return _build_blast_result(
749
+ hits,
750
+ source="ncbi",
751
+ database=database,
752
+ program=program,
753
+ seq_type=seq_type,
754
+ query_accession=query_accession,
755
+ query_length=parsed.get("query_length", 0),
756
+ )
757
+ ncbi_error = parsed["error"]
758
+ else:
759
+ ncbi_error = results["error"]
760
+
761
+ logger.warning("EBI BLAST unavailable and NCBI failed (%s)", ncbi_error)
762
+ return {"error": ncbi_error or "BLAST failed via EBI and NCBI", "count": 0, "hits": []}
763
+
764
+
765
+ async def _run_uniprot(top_hit: dict, query_sequence: str | None = None, try_sequence: bool = True) -> dict:
766
+ accession = (top_hit.get("accession") or "").strip()
767
+ if not accession:
768
+ return {"error": "No accession"}
769
+ desc = (top_hit.get("description") or "").strip()
770
+ organism = (top_hit.get("organism") or "").strip()
771
+
772
+ # Resolve ANY BLAST-hit accession (RefSeq, GenBank, PDB, Ensembl, ...) to a
773
+ # UniProt accession via the strategy ladder in identifier_resolution:
774
+ # direct -> xref search (fast) -> name search -> EBI sequence BLAST -> idmapping
775
+ # Exhausting all tiers returns an explicit unresolved/de_novo result (§1).
776
+ try:
777
+ resolved = await resolve_to_uniprot(
778
+ accession=accession,
779
+ sequence=query_sequence,
780
+ description=desc,
781
+ organism=organism or None,
782
+ try_sequence=try_sequence,
783
+ )
784
+ except Exception as e:
785
+ logger.warning("Identifier resolution failed for %s: %s", accession, e)
786
+ resolved = None
787
+
788
+ if not resolved or resolved.get("status") != "resolved":
789
+ logger.info("No UniProt mapping for %s, using BLAST data only", accession)
790
+ return {
791
+ "accession": accession,
792
+ "full_name": top_hit.get("description", ""),
793
+ "organism": top_hit.get("organism", ""),
794
+ "gene_names": [],
795
+ "functions": [],
796
+ "keywords": [],
797
+ "subcellular_locations": [],
798
+ "pdb_ids": [],
799
+ "go_terms": [],
800
+ "sequence": "",
801
+ "sequence_length": 0,
802
+ "features": [],
803
+ "resolution": {
804
+ "uniprot_accession": None,
805
+ "method": "de_novo",
806
+ "original_accession": accession,
807
+ },
808
+ "resolved_uniprot": False,
809
+ # BLAST-derived similarity is homolog-grade certainty at best (§1.1)
810
+ "confidence": "homolog",
811
+ "_note": f"No UniProt mapping found for {accession} — showing BLAST-derived data only",
812
+ }
813
+
814
+ uniprot_acc = resolved["accession"]
815
+ method = resolved["method"]
816
+
817
+ tool = UniprotTool()
818
+ result = await tool.run({"accession": uniprot_acc})
819
+ if "error" in result:
820
+ return {"error": result["error"]}
821
+
822
+ return {
823
+ "accession": result.get("accession", uniprot_acc),
824
+ "full_name": result.get("full_name", ""),
825
+ "organism": result.get("organism", ""),
826
+ "gene_names": result.get("gene_names", []),
827
+ "functions": result.get("functions", []),
828
+ "keywords": result.get("keywords", []),
829
+ "subcellular_locations": result.get("subcellular_locations", []),
830
+ "pdb_ids": result.get("pdb_ids", []),
831
+ "go_terms": result.get("go_terms", []),
832
+ "sequence": result.get("sequence", ""),
833
+ "sequence_length": result.get("sequence_length", 0),
834
+ "features": [
835
+ f for f in (result.get("features", []) or [])
836
+ if f.get("type") in ("ACTIVE_SITE", "BINDING", "MUTAGENESIS", "SITE", "MOD_RES")
837
+ ],
838
+ "resolution": {
839
+ "uniprot_accession": uniprot_acc,
840
+ "method": method,
841
+ "original_accession": accession,
842
+ },
843
+ "resolved_uniprot": True,
844
+ "confidence": resolved.get("confidence", "identified"),
845
+ }
846
+
847
+
848
+ async def _run_msa(query_sequence: str, blast_hits: list, alignment_mode: str = "global") -> dict:
849
+ sequences = [("query", query_sequence)]
850
+
851
+ for hit in blast_hits[:5]:
852
+ acc = hit.get("accession", "")
853
+ hit_seq = hit.get("hit_alignment", "")
854
+
855
+ if acc:
856
+ source = detect_source_from_accession(acc)
857
+ mapped_acc = acc
858
+ if source == "ncbi":
859
+ mapped = await map_refseq_to_uniprot(acc)
860
+ if mapped:
861
+ mapped_acc = mapped
862
+ try:
863
+ tool = UniprotTool()
864
+ ud = await tool.run({"accession": mapped_acc})
865
+ if "error" not in ud and ud.get("sequence"):
866
+ clean_seq = "".join(c for c in ud["sequence"] if c.isalpha()).upper()
867
+ if len(clean_seq) > 10:
868
+ sequences.append((acc, clean_seq))
869
+ continue
870
+ except Exception:
871
+ pass
872
+
873
+ if hit_seq:
874
+ clean = "".join(c for c in hit_seq if c.isalpha()).upper()
875
+ if len(clean) > 10:
876
+ sequences.append((f"{acc}_aln", clean))
877
+
878
+ if len(sequences) < 2:
879
+ return {"error": "Not enough sequences for MSA", "aln_fasta": None, "phylotree": None}
880
+
881
+ fasta_lines = []
882
+ for sid, sseq in sequences:
883
+ fasta_lines.append(f">{sid}")
884
+ for i in range(0, len(sseq), 80):
885
+ fasta_lines.append(sseq[i:i + 80])
886
+ fasta_str = "\n".join(fasta_lines)
887
+
888
+ try:
889
+ email = settings.NCBI_EMAIL or "bioflow@example.com"
890
+ seq_type = detect_sequence_type(query_sequence) or "protein"
891
+ stype = "protein" if seq_type == "protein" else "dna"
892
+ try:
893
+ # Try local MAFFT first (fast, no network dependency)
894
+ from app.tools.mafft_local import run_local_mafft
895
+ local_result = await asyncio.get_running_loop().run_in_executor(
896
+ None, run_local_mafft, fasta_str, "auto", 1, 300,
897
+ )
898
+ if local_result and local_result.get("aln_fasta"):
899
+ method = "mafft-local"
900
+ aln_fasta = local_result["aln_fasta"]
901
+ phylotree = ""
902
+ else:
903
+ raise ValueError("local MAFFT unavailable or returned empty")
904
+ except Exception:
905
+ try:
906
+ result = await run_ebi_msa_best_effort(
907
+ sequence=fasta_str,
908
+ stype=stype,
909
+ email=email,
910
+ )
911
+ method = result["method"]
912
+ aln_fasta = result["aln_fasta"]
913
+ phylotree = result["phylotree"]
914
+ except Exception as e:
915
+ logger.warning("EBI MSA unavailable (%s) — using in-process fallback", e)
916
+ from app.tools.msa_fallback import progressive_msa
917
+ fallback = await asyncio.get_running_loop().run_in_executor(
918
+ None, progressive_msa, sequences, stype
919
+ )
920
+ aln_fasta, phylotree = fallback
921
+ method = "in-process fallback"
922
+
923
+ payload = {
924
+ "aln_fasta": aln_fasta,
925
+ "phylotree": phylotree,
926
+ "sequence_count": len(sequences),
927
+ "alignment_mode": alignment_mode,
928
+ "method": method,
929
+ "_fallback": method != "clustalo",
930
+ }
931
+
932
+ # Local mode: refine query vs the best non-query sequence with an
933
+ # in-process Smith-Waterman alignment so the wizard can show which
934
+ # region actually matches (EBI MSA itself is always global).
935
+ if alignment_mode == "local" and len(sequences) > 1:
936
+ subject_id, subject_seq = sequences[1]
937
+ try:
938
+ payload["pairwise"] = pairwise_align(query_sequence, subject_seq, mode="local")
939
+ payload["pairwise_subject"] = subject_id
940
+ except Exception as e:
941
+ logger.warning("Local pairwise refinement failed: %s", e)
942
+
943
+ return payload
944
+
945
+ except Exception as e:
946
+ return {"error": str(e), "aln_fasta": None, "phylotree": None}
947
+
948
+
949
+ async def _run_domains(accession: str) -> dict:
950
+ """Run domain analysis using the shared tool module (eliminates code duplication)."""
951
+ try:
952
+ from app.tools.domain_analysis import fetch_interpro_domains
953
+ return await fetch_interpro_domains(accession)
954
+ except Exception as e:
955
+ return {"error": str(e), "uniprot_accession": accession, "sequence_length": 0, "domains": []}
956
+
957
+
958
+ async def _run_denovo_steps(job_id: str, sequence: str, steps: list[str], _mark, context: dict) -> None:
959
+ """Tier-6 branch (techspec.md §1): characterize from sequence alone.
960
+
961
+ Runs when BLAST finds no homolog at all. Composition/function hints land
962
+ in the uniprot slot (the annotation slot); MSA/phylo/pathways are marked
963
+ explicitly unavailable rather than left empty.
964
+ """
965
+ import asyncio as _asyncio
966
+ from app.services.de_novo import (
967
+ composition_stats, esmfold_structure, function_hints, interpro_sequence_search,
968
+ )
969
+
970
+ unavailable = "Unavailable for de novo sequences — no identified homolog"
971
+
972
+ async def _do_annotation():
973
+ bundle = {"_de_novo": True}
974
+ try:
975
+ bundle["composition"] = composition_stats(sequence)
976
+ except Exception as e:
977
+ bundle["composition"] = {"error": str(e)}
978
+ try:
979
+ bundle["function_hints"] = function_hints(sequence)
980
+ except Exception as e:
981
+ bundle["function_hints"] = {"error": str(e)}
982
+ return bundle
983
+
984
+ async def _do_domains():
985
+ try:
986
+ return await interpro_sequence_search(sequence)
987
+ except Exception as e:
988
+ logger.warning("De novo InterProScan failed: %s", e)
989
+ return {"error": str(e), "domains": [], "sequence_length": 0}
990
+
991
+ async def _do_structure():
992
+ try:
993
+ return await esmfold_structure(sequence)
994
+ except Exception as e:
995
+ logger.warning("De novo ESMFold failed: %s", e)
996
+ return {"structure_available": False, "source": "esmfold", "message": str(e)}
997
+
998
+ fan = []
999
+ names = []
1000
+ if "uniprot" in steps:
1001
+ fan.append(_do_annotation())
1002
+ names.append("uniprot")
1003
+ if "domains" in steps:
1004
+ fan.append(_do_domains())
1005
+ names.append("domains")
1006
+ if "alphafold" in steps:
1007
+ fan.append(_do_structure())
1008
+ names.append("alphafold")
1009
+
1010
+ for name in names:
1011
+ _mark(name, "running", progress=10)
1012
+ results = await _asyncio.gather(*fan, return_exceptions=True)
1013
+
1014
+ for name, res in zip(names, results):
1015
+ if isinstance(res, Exception):
1016
+ _mark(name, "failed", error=str(res)[:500])
1017
+ continue
1018
+ ok = (
1019
+ res.get("domains") is not None if name == "domains"
1020
+ else res.get("structure_available") is True if name == "alphafold"
1021
+ else True
1022
+ )
1023
+ _mark(name, "complete" if ok else "failed", progress=100,
1024
+ data=res, error=None if ok else "Prediction returned nothing usable")
1025
+ context[name] = res
1026
+
1027
+ # Annotation-database features have no de novo substitute — say so.
1028
+ for name in ("msa", "phylo", "pathway_enrichment"):
1029
+ if name in steps:
1030
+ _mark(name, "failed", error=unavailable)
1031
+
1032
+
1033
+ async def _run_domains_or_denovo(sequence: str, accession: str | None, resolved_uniprot: bool) -> dict:
1034
+ """Accession lookup when resolved; InterProScan sequence-search otherwise."""
1035
+ if accession and resolved_uniprot:
1036
+ return await _run_domains(accession)
1037
+ from app.services.de_novo import interpro_sequence_search
1038
+ try:
1039
+ result = await interpro_sequence_search(sequence)
1040
+ result["confidence"] = "de_novo"
1041
+ return result
1042
+ except Exception as e:
1043
+ return {"error": str(e), "domains": [], "sequence_length": 0}
1044
+
1045
+
1046
+ async def _run_alphafold_or_esmfold(
1047
+ context: dict, sequence: str, accession: str | None, resolved_uniprot: bool,
1048
+ ) -> dict:
1049
+ """AlphaFold DB lookup when resolved; ESMFold ab initio otherwise."""
1050
+ if accession and resolved_uniprot:
1051
+ result = await _run_alphafold(context)
1052
+ return result or {}
1053
+ from app.services.de_novo import esmfold_structure
1054
+ try:
1055
+ return await esmfold_structure(sequence)
1056
+ except Exception as e:
1057
+ return {"structure_available": False, "source": "esmfold", "message": str(e)}
1058
+
1059
+
1060
+ async def _run_pathway_enrichment(context: dict) -> dict | None:
1061
+ gene_names = []
1062
+ uniprot = context.get("uniprot", {})
1063
+ if isinstance(uniprot, dict):
1064
+ gene_names = uniprot.get("gene_names", [])[:20] if isinstance(uniprot.get("gene_names"), list) else []
1065
+ if not gene_names:
1066
+ blast_data = context.get("blast", {})
1067
+ if isinstance(blast_data, dict):
1068
+ for hit in (blast_data.get("hits") or [])[:10]:
1069
+ words = (hit.get("description", "") or "").replace("(", " ").replace(")", " ").split()
1070
+ for w in words:
1071
+ if w.isupper() and len(w) >= 2 and not w.startswith("OS="):
1072
+ gene_names.append(w)
1073
+ break
1074
+ if not gene_names:
1075
+ return None
1076
+ try:
1077
+ from app.services.pathway_enrichment import run_enrichment
1078
+ result = await run_enrichment(gene_names)
1079
+ return result
1080
+ except Exception as e:
1081
+ logger.warning(f"Pathway enrichment failed: {e}")
1082
+ return None
1083
+
1084
+
1085
+ async def _run_alphafold(context: dict) -> dict | None:
1086
+ uniprot_data = context.get("uniprot", {})
1087
+ accession = uniprot_data.get("accession") if isinstance(uniprot_data, dict) else None
1088
+ if not accession:
1089
+ return None
1090
+ try:
1091
+ from app.tools.alphafold import AlphaFoldTool
1092
+ result = await AlphaFoldTool().run({"uniprot_accession": accession})
1093
+ return result
1094
+ except Exception as e:
1095
+ logger.warning(f"AlphaFold fetch failed for {accession}: {e}")
1096
+ return {"structure_available": False, "message": str(e)}
1097
+
1098
+
1099
+ async def _run_interpret(context: dict) -> dict:
1100
+ from app.ai.interpreter import interpret_text
1101
+ prompt_context = {
1102
+ "blast": context.get("blast", {}),
1103
+ "uniprot": context.get("uniprot", {}),
1104
+ "alphafold": context.get("alphafold", {}),
1105
+ "pathway_enrichment": context.get("pathway_enrichment", {}),
1106
+ "query_confidence": context.get("query", {}).get("confidence", "identified"),
1107
+ }
1108
+ return await interpret_text("protein_analysis", prompt_context)
migrations/011_ngs_production_runs.sql ADDED
@@ -0,0 +1,26 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ -- Durable ownership and execution state for real nf-core/sarek runs.
2
+ CREATE TABLE IF NOT EXISTS ngs_production_runs (
3
+ run_id uuid PRIMARY KEY,
4
+ user_id uuid NOT NULL REFERENCES auth.users(id) ON DELETE CASCADE,
5
+ state text NOT NULL CHECK (state IN ('SUBMITTED','PENDING','RUNNING','SUCCEEDED','FAILED','UNKNOWN')),
6
+ executor text NOT NULL CHECK (executor IN ('local','slurm','awsbatch')),
7
+ executor_job_id text NOT NULL,
8
+ workflow text NOT NULL,
9
+ revision text NOT NULL,
10
+ outdir text NOT NULL,
11
+ submitted_at timestamptz NOT NULL,
12
+ updated_at timestamptz NOT NULL,
13
+ exit_code integer,
14
+ message text
15
+ );
16
+
17
+ CREATE INDEX IF NOT EXISTS idx_ngs_production_runs_user_updated
18
+ ON ngs_production_runs(user_id, updated_at DESC);
19
+
20
+ ALTER TABLE ngs_production_runs ENABLE ROW LEVEL SECURITY;
21
+
22
+ DO $$ BEGIN
23
+ CREATE POLICY "Users can view own production NGS runs"
24
+ ON ngs_production_runs FOR SELECT USING (auth.uid() = user_id);
25
+ EXCEPTION WHEN duplicate_object THEN NULL;
26
+ END $$;
requirements.txt CHANGED
@@ -21,3 +21,4 @@ pymol-open-source-whl>=3.2 # unofficial wheel of open-source PyMOL; provides th
21
  # dev / testing
22
  pytest>=7.0
23
  httpx>=0.24
 
 
21
  # dev / testing
22
  pytest>=7.0
23
  httpx>=0.24
24
+ boto3>=1.35,<2
tests/test_ngs_platform.py CHANGED
@@ -545,6 +545,15 @@ def test_reference_validation_resolves_grch38():
545
  assert out["summary"]["reference"]["build"] == "GRCh38"
546
 
547
 
 
 
 
 
 
 
 
 
 
548
  def test_reference_validation_mismatch_stops():
549
  # GRCh38 data against a GRCh37 annotation is a provenance error -> STOP.
550
  out = run_reference_validation({"reference": "grch38", "annotation_build": "GRCh37"})
@@ -623,7 +632,9 @@ def test_choose_aligner_by_assay():
623
  def test_alignment_run_and_mapping_ok():
624
  ref, reads = _make_ref_and_reads(ref_len=600, n_reads=20)
625
  out = run_alignment(reads[:20], ref, assay="WGS", read_length=30)
626
- assert out["summary"]["meta"]["aligner"] == "bwa-mem2"
 
 
627
  assert out["summary"]["decision"] != "STOP"
628
  assert out["summary"]["meta"]["mapped"] >= 18
629
  assert len(out["records"]) == 20
 
545
  assert out["summary"]["reference"]["build"] == "GRCh38"
546
 
547
 
548
+ def test_synthetic_reference_never_claims_grch38_coordinates():
549
+ out = run_reference_validation({"reference": "grch38", "synthetic_reference": True})
550
+ reference = out["summary"]["reference"]
551
+ assert out["summary"]["status"] == "PASS"
552
+ assert reference["id"] == "synthetic-positive-control"
553
+ assert reference["build"] == "NOT_APPLICABLE"
554
+ assert reference["requested_template"] == "grch38"
555
+
556
+
557
  def test_reference_validation_mismatch_stops():
558
  # GRCh38 data against a GRCh37 annotation is a provenance error -> STOP.
559
  out = run_reference_validation({"reference": "grch38", "annotation_build": "GRCh37"})
 
632
  def test_alignment_run_and_mapping_ok():
633
  ref, reads = _make_ref_and_reads(ref_len=600, n_reads=20)
634
  out = run_alignment(reads[:20], ref, assay="WGS", read_length=30)
635
+ assert out["summary"]["meta"]["recommended_production_aligner"] == "bwa-mem2"
636
+ assert out["summary"]["meta"]["executed_implementation"] == "bionexus-seed-aligner"
637
+ assert out["summary"]["meta"]["executor"] == "python-surrogate"
638
  assert out["summary"]["decision"] != "STOP"
639
  assert out["summary"]["meta"]["mapped"] >= 18
640
  assert len(out["records"]) == 20
tests/test_ngs_production_execution.py ADDED
@@ -0,0 +1,45 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from pathlib import Path
2
+
3
+ from app.config import settings
4
+ from app.ngs import execution
5
+
6
+
7
+ def test_capabilities_never_offer_preview_fallback(monkeypatch):
8
+ monkeypatch.setattr(settings, "NGS_LOCAL_EXECUTION_ENABLED", False)
9
+ monkeypatch.setattr(settings, "NGS_SLURM_EXECUTION_ENABLED", False)
10
+ monkeypatch.setattr(settings, "NGS_AWS_BATCH_EXECUTION_ENABLED", False)
11
+ result = execution.executor_capabilities()
12
+ assert result["fallback"] is None
13
+ assert all(not item["available"] for item in result["executors"].values())
14
+ assert "never fall back" in result["note"]
15
+
16
+
17
+ def test_disabled_executor_fails_closed(monkeypatch):
18
+ monkeypatch.setattr(settings, "NGS_LOCAL_EXECUTION_ENABLED", False)
19
+ try:
20
+ execution.submit_run("local", ["nextflow", "run", "nf-core/sarek"], "/results", "user-1")
21
+ except RuntimeError as exc:
22
+ assert "executor unavailable" in str(exc)
23
+ else:
24
+ raise AssertionError("disabled executor accepted a production run")
25
+
26
+
27
+ def test_owned_run_record_round_trip(tmp_path, monkeypatch):
28
+ monkeypatch.setattr(settings, "NGS_RUN_ROOT", str(tmp_path))
29
+ monkeypatch.setattr("app.services.supabase.get_supabase", lambda: (_ for _ in ()).throw(RuntimeError("offline")))
30
+ record = {
31
+ "run_id": "11111111-1111-1111-1111-111111111111",
32
+ "state": "SUBMITTED", "executor": "local", "executor_job_id": "123",
33
+ "workflow": "nf-core/sarek", "revision": "3.10.0", "outdir": "/results",
34
+ "submitted_at": "2026-09-02T00:00:00+00:00", "updated_at": "2026-09-02T00:00:00+00:00",
35
+ "exit_code": None, "message": None, "user_id": "user-1",
36
+ }
37
+ execution._persist_record(record)
38
+ assert Path(tmp_path, f"{record['run_id']}.json").is_file()
39
+ assert execution._load_record(record["run_id"], "user-1")["workflow"] == "nf-core/sarek"
40
+ try:
41
+ execution._load_record(record["run_id"], "user-2")
42
+ except FileNotFoundError:
43
+ pass
44
+ else:
45
+ raise AssertionError("run record was visible to another user")
tests/test_ngs_v2_router.py CHANGED
@@ -1,4 +1,6 @@
1
  import gzip
 
 
2
  import os
3
  import random
4
 
@@ -7,6 +9,9 @@ from fastapi import FastAPI
7
  from fastapi.testclient import TestClient
8
 
9
  from app.routers import ngs_v2
 
 
 
10
 
11
 
12
  @pytest.fixture(scope="module")
@@ -40,6 +45,167 @@ def test_stages_lists_contracts(client):
40
  assert len(steps) == 21
41
 
42
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
43
  def test_detect_returns_evidence(client, tmp_path):
44
  path, _ = _write_fastq(tmp_path, "HUM0001_R1.fastq.gz", 40)
45
  r = client.post("/api/ngs/v2/detect", json={"file_paths": [path]})
@@ -90,6 +256,34 @@ def test_analyze_runs_full_dag_through_final_gate(client, tmp_path):
90
  assert len(provenance["inputs"]) == 2
91
  assert all(item["checksum"]["algorithm"] == "md5" for item in provenance["inputs"])
92
  assert len(provenance["tools"]) == 21
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
93
 
94
 
95
  def test_analyze_emits_igv_tracks(client, tmp_path):
 
1
  import gzip
2
+ import hashlib
3
+ import hmac
4
  import os
5
  import random
6
 
 
9
  from fastapi.testclient import TestClient
10
 
11
  from app.routers import ngs_v2
12
+ from app.models.responses import NgsClinicalEvidenceRequest
13
+ from app.ngs.production import clinical_signature_payload
14
+ from app.config import settings
15
 
16
 
17
  @pytest.fixture(scope="module")
 
45
  assert len(steps) == 21
46
 
47
 
48
+ def test_portable_benchmark_reports_scoped_execution_parity(client):
49
+ response = client.get("/api/ngs/v2/benchmarks/portable")
50
+ assert response.status_code == 200
51
+ report = response.json()
52
+ assert report["status"] == "PASS"
53
+ assert report["workflow_output_parity"] is True
54
+ assert report["expected_call"] == {
55
+ "chrom": "chrTiny", "pos": 50, "ref": "C", "alt": "G",
56
+ "genotype": "0/1", "depth": 20, "allelic_depth": "10,10",
57
+ }
58
+ assert len({row["normalized_sha256"] for row in report["reports"]}) == 1
59
+ assert all(row["f1"] == 1.0 for row in report["reports"])
60
+ galaxy = next(row for row in report["reports"] if row["orchestrator"].startswith("Galaxy"))
61
+ assert galaxy["execution"] == "EXECUTED_WITHOUT_GALAXY_SERVER"
62
+
63
+
64
+ def test_production_capabilities_are_explicit_and_have_no_preview_fallback(client):
65
+ response = client.get("/api/ngs/v2/production/capabilities")
66
+ assert response.status_code == 200
67
+ body = response.json()
68
+ assert set(body["executors"]) == {"local", "slurm", "awsbatch"}
69
+ assert body["fallback"] is None
70
+ assert body["workflow"] == {"name": "nf-core/sarek", "revision": "3.10.0"}
71
+
72
+
73
+ def test_clean_demo_uses_explicit_not_evaluated_states(client):
74
+ response = client.post("/api/ngs/v2/analyze", json={"demo_profile": "wgs-clean"})
75
+ assert response.status_code == 200
76
+ body = response.json()
77
+ stages = {stage["step"]: stage for stage in body["pipeline"]["stages"]}
78
+
79
+ assert body["requested"]["reference"] == "synthetic-positive-control"
80
+ assert body["requested"]["reference_template_requested"] == "grch38"
81
+ assert body["detection"]["sample_type"] == "synthetic-positive-control"
82
+ reference = body["pipeline"]["provenance"]["reference"]
83
+ assert reference["id"] == "synthetic-positive-control"
84
+ assert reference["reference_kind"] == "SYNTHETIC"
85
+ assert reference["build"] == "NOT_APPLICABLE"
86
+ assert reference["artifacts"] == [{"name": "In-memory synthetic reference", "present": True}]
87
+
88
+ assert stages["coverage"]["qc"]["status"] == "PASS"
89
+ assert stages["variant_filter"]["qc"]["status"] == "PASS"
90
+ assert stages["variant_filter"]["data"]["n_final"] == 0
91
+ assert stages["annotation"]["qc"]["status"] == "PASS"
92
+ assert stages["annotation"]["data"]["n_annotated"] == 0
93
+ assert stages["annotation"]["data"]["status"] == "NOT_APPLICABLE"
94
+ assert stages["contamination"]["data"]["status"] == "NOT_EVALUATED"
95
+ assert stages["contamination"]["qc"]["metrics"][0]["value"] is None
96
+ assert stages["identity"]["data"]["status"] == "NOT_EVALUATED"
97
+ assert len(stages["identity"]["qc"]["metrics"]) == 1
98
+ assert stages["final_gate"]["qc"]["metrics"][0]["value"] == "ANALYSIS_READY_WITH_WARNINGS"
99
+
100
+ warnings = body["pipeline"]["warnings"]
101
+ assert len(warnings) == 2
102
+ assert all("Not evaluated:" in warning for warning in warnings)
103
+ assert not any("(None)" in warning for warning in warnings)
104
+ assert not any(name in " ".join(warnings) for name in (
105
+ "contam_checked", "identity_checked", "candidate_frac", "gate_ok",
106
+ ))
107
+ assert all(tool["implementation"] for tool in body["pipeline"]["provenance"]["tools"])
108
+
109
+
110
+ def test_production_plan_pins_sarek_and_emits_auditable_argv(client):
111
+ response = client.post("/api/ngs/v2/production/plan", json={
112
+ "assay": "WGS",
113
+ "sample_model": "trio",
114
+ "input_type": "FASTQ",
115
+ "samplesheet_path": "/staged/family.csv",
116
+ "outdir": "/results/family",
117
+ "genome": "GRCh38",
118
+ "execution_profile": "docker",
119
+ "caller": "haplotypecaller",
120
+ "clinical_intent": True,
121
+ })
122
+ assert response.status_code == 200
123
+ plan = response.json()
124
+ assert plan["ready_to_launch"] is True
125
+ assert plan["state"] == "PLANNED"
126
+ assert plan["workflow"]["name"] == "nf-core/sarek"
127
+ assert plan["workflow"]["revision"] == "3.10.0"
128
+ assert plan["command_argv"][:5] == ["nextflow", "run", "nf-core/sarek", "-r", "3.10.0"]
129
+ assert ["--genome", "GATK.GRCh38"] == plan["command_argv"][plan["command_argv"].index("--genome"):plan["command_argv"].index("--genome") + 2]
130
+ assert "--joint_germline" in plan["command_argv"]
131
+ assert plan["clinical_boundary"]["current_status"] == "NOT_CLINICALLY_RELEASABLE"
132
+ assert {item["id"] for item in plan["required_artifacts"]} >= {"execution", "multiqc", "alignment", "small_variants", "provenance"}
133
+
134
+
135
+ def test_production_wes_and_remote_executor_fail_closed(client):
136
+ response = client.post("/api/ngs/v2/production/plan", json={
137
+ "assay": "WES",
138
+ "samplesheet_path": "/staged/case.csv",
139
+ "outdir": "/results/case",
140
+ "execution_profile": "slurm",
141
+ })
142
+ assert response.status_code == 200
143
+ plan = response.json()
144
+ assert plan["state"] == "BLOCKED"
145
+ assert plan["ready_to_launch"] is False
146
+ assert any("target BED" in item for item in plan["blockers"])
147
+ assert any("custom configuration" in item for item in plan["blockers"])
148
+
149
+
150
+ def test_clinical_gate_rejects_claims_without_evidence(client):
151
+ response = client.post("/api/ngs/v2/clinical/evaluate", json={})
152
+ assert response.status_code == 200
153
+ result = response.json()
154
+ assert result["status"] == "NOT_CLINICALLY_RELEASABLE"
155
+ assert result["clinically_validated"] is False
156
+ assert result["missing_or_failed"]
157
+ signature_gate = next(gate for gate in result["gates"] if gate["id"] == "evidence_signature")
158
+ assert signature_gate["status"] == "FAIL"
159
+
160
+
161
+ def test_clinical_gate_passes_only_complete_signed_external_evidence(client, monkeypatch):
162
+ sha = "a" * 64
163
+ payload = {
164
+ "evidence_bundle_sha256": "d" * 64,
165
+ "assay_validation_id": "VAL-WGS-2026-04",
166
+ "workflow_status": "COMPLETED",
167
+ "sarek_revision": "3.10.0",
168
+ "reference_build": "GRCh38",
169
+ "reference_manifest_sha256": sha,
170
+ "samplesheet_sha256": "b" * 64,
171
+ "container_digests_complete": True,
172
+ "complete_input_processed": True,
173
+ "required_artifacts_present": True,
174
+ "qc_pass": True,
175
+ "sample_identity_pass": True,
176
+ "contamination_pass": True,
177
+ "sex_ploidy_reviewed": True,
178
+ "truthset_name": "GIAB HG002 v4.2.1",
179
+ "benchmark_protocol_id": "SOP-BENCH-07",
180
+ "benchmark_acceptance_pass": True,
181
+ "confident_regions_sha256": "c" * 64,
182
+ "same_sample_reference_regions": True,
183
+ "snv_precision": 0.999,
184
+ "snv_recall": 0.998,
185
+ "indel_precision": 0.995,
186
+ "indel_recall": 0.994,
187
+ "human_reviewed": True,
188
+ "reviewer_id": "reviewer-17",
189
+ "release_signature_id": "sig-2026-09-01-001",
190
+ }
191
+ signing_key = "test-only-clinical-attestation-key"
192
+ monkeypatch.setattr(settings, "NGS_CLINICAL_EVIDENCE_HMAC_KEY", signing_key)
193
+ evidence = NgsClinicalEvidenceRequest(**payload)
194
+ payload["evidence_signature"] = hmac.new(
195
+ signing_key.encode("utf-8"), clinical_signature_payload(evidence), hashlib.sha256,
196
+ ).hexdigest()
197
+ response = client.post("/api/ngs/v2/clinical/evaluate", json=payload)
198
+ assert response.status_code == 200
199
+ result = response.json()
200
+ assert result["status"] == "SOFTWARE_GATE_PASSED"
201
+ assert result["clinically_validated"] is False
202
+ assert result["missing_or_failed"] == []
203
+
204
+ payload["truthset_name"] = "synthetic positive control"
205
+ rejected = client.post("/api/ngs/v2/clinical/evaluate", json=payload).json()
206
+ assert rejected["status"] == "NOT_CLINICALLY_RELEASABLE"
207
+
208
+
209
  def test_detect_returns_evidence(client, tmp_path):
210
  path, _ = _write_fastq(tmp_path, "HUM0001_R1.fastq.gz", 40)
211
  r = client.post("/api/ngs/v2/detect", json={"file_paths": [path]})
 
256
  assert len(provenance["inputs"]) == 2
257
  assert all(item["checksum"]["algorithm"] == "md5" for item in provenance["inputs"])
258
  assert len(provenance["tools"]) == 21
259
+ assert all("implementation" in item for item in provenance["tools"])
260
+ assert all("evidence_level" in stage for stage in stages)
261
+
262
+ validation = body["pipeline"]["validation"]
263
+ assert validation["claim"] == "NO_ACCURACY_CLAIM"
264
+ assert validation["same_or_better_supported"] is False
265
+ assert len(validation["comparisons"]) >= 3
266
+ assert all(item["status"] != "EVALUATED" for item in validation["comparisons"])
267
+ assert all(item["metrics"] is None for item in validation["comparisons"])
268
+ assert validation["analysis_grade"] == "EXPLORATORY_PREVIEW"
269
+ assert validation["research_ready"] is False
270
+ assert all(item["status"] == "MISSING" for item in validation["production_requirements"])
271
+ assert body["requested"]["all_records_processed"] is True
272
+
273
+
274
+ def test_analyze_discloses_fastq_sampling(client, tmp_path):
275
+ r1, _ = _write_fastq(tmp_path, "HG002_R1.fastq.gz", 2001, seed=21)
276
+ r2, _ = _write_fastq(tmp_path, "HG002_R2.fastq.gz", 2001, seed=22)
277
+ response = client.post("/api/ngs/v2/analyze", json={
278
+ "file_paths": [r1, r2], "reference": "grch38", "assay": "WGS",
279
+ "synthetic_reference": True,
280
+ })
281
+ assert response.status_code == 200
282
+ body = response.json()
283
+ assert body["requested"]["record_cap_per_file"] == 2000
284
+ assert body["requested"]["all_records_processed"] is False
285
+ assert sorted(body["requested"]["truncated_files"]) == ["HG002_R1.fastq.gz", "HG002_R2.fastq.gz"]
286
+ assert body["pipeline"]["validation"]["input_sampling"]["mode"] == "SAMPLED_PREVIEW"
287
 
288
 
289
  def test_analyze_emits_igv_tracks(client, tmp_path):