Samad14 commited on
Commit
d03beec
·
verified ·
1 Parent(s): 89ce8b4

Upload folder using huggingface_hub

Browse files
app/rnaseq/deseq2_analysis.R CHANGED
@@ -54,6 +54,11 @@ if (nzchar(covariates_raw)) {
54
  if (any(!vapply(covariates, safe_name, logical(1)))) stop("covariate name has unsupported characters")
55
  }
56
 
 
 
 
 
 
57
  txi <- NULL
58
  if (input_kind == "salmon") {
59
  suppressPackageStartupMessages(library(tximport))
@@ -85,6 +90,8 @@ if (input_kind == "salmon") {
85
  }
86
  if (anyDuplicated(colnames(counts)) > 0) stop("count matrix contains duplicate sample names")
87
 
 
 
88
  meta <- read.delim(metadata_path, check.names = FALSE, stringsAsFactors = FALSE)
89
  if (!("sample" %in% names(meta))) stop("metadata must contain a 'sample' column")
90
  if (!(condition_col %in% names(meta))) stop(paste0("metadata is missing condition column '", condition_col, "'"))
@@ -241,17 +248,23 @@ design_audit <- list(
241
  if (length(design_warnings) == 0) design_audit$status <- "PASS"
242
  write(toJSON(design_audit, auto_unbox = TRUE, pretty = TRUE, null = "null", digits = 10), file.path(outdir, "design_audit.json"))
243
 
 
 
244
  dds <- if (is.null(txi)) {
245
  DESeqDataSetFromMatrix(countData = counts, colData = meta, design = design_formula)
246
  } else {
247
  DESeqDataSetFromTximport(txi, colData = meta, design = design_formula)
248
  }
 
 
249
  genes_input <- nrow(dds)
250
  keep <- rowSums(counts(dds) >= min_count) >= min_samples
251
  dds <- dds[keep, ]
252
  genes_kept <- nrow(dds)
253
  if (genes_kept < 2) stop("pre-filtering retained fewer than two genes")
254
 
 
 
255
  dds <- estimateSizeFactors(dds)
256
  normalized <- counts(dds, normalized = TRUE)
257
  size_factors <- sizeFactors(dds)
@@ -275,6 +288,8 @@ size_factor_df <- data.frame(sample = names(size_factors), size_factor = as.nume
275
  write.table(size_factor_df, file.path(outdir, "size_factors.tsv"), sep = "\t", quote = FALSE, row.names = FALSE)
276
  write.table(data.frame(gene = rownames(normalized), normalized, check.names = FALSE), file.path(outdir, "normalized_counts.tsv"), sep = "\t", quote = FALSE, row.names = FALSE)
277
 
 
 
278
  # vst() is fast for ordinary whole-transcriptome matrices but requires enough rows
279
  # for its trend subsampling. varianceStabilizingTransformation() is the exact
280
  # DESeq2 transform and is robust for small deterministic validation subsets.
@@ -325,6 +340,8 @@ if (stage == "qc") {
325
  quit(save = "no", status = 0)
326
  }
327
 
 
 
328
  dds <- DESeq(dds)
329
  res <- results(dds, contrast = c(condition_col, test_level, reference_level), alpha = alpha)
330
  res_df <- as.data.frame(res)
@@ -358,6 +375,8 @@ write.table(res_df, file.path(outdir, "deseq2_all_results.tsv"), sep = "\t", quo
358
  write.table(deg_df, file.path(outdir, "deseq2_significant.tsv"), sep = "\t", quote = FALSE, row.names = FALSE, na = "")
359
 
360
 
 
 
361
  pca$.condition <- pca[[condition_col]]
362
  pca_plot <- ggplot(pca, aes(x = PC1, y = PC2, label = name, shape = .condition)) +
363
  geom_point(size = 3.2) +
@@ -399,6 +418,8 @@ volcano_plot <- ggplot(volcano_df, aes(x = log2FoldChange, y = minus_log10_padj,
399
  theme_minimal(base_size = 11) + theme(legend.position = "bottom")
400
  plot_svg_pdf_png("volcano", 8, 6, function() print(volcano_plot))
401
 
 
 
402
  heatmap_written <- FALSE
403
  heatmap_basis <- "none"
404
  heatmap_gene_count <- 0L
@@ -445,6 +466,8 @@ if (length(selected_genes) >= 2) {
445
  heatmap_written <- TRUE
446
  }
447
 
 
 
448
  script_arg <- grep("^--file=", commandArgs(), value = TRUE)[[1]]
449
  source(file.path(dirname(sub("^--file=", "", script_arg)), "go_enrichment.R"))
450
  enrichment <- tryCatch(
@@ -458,6 +481,7 @@ enrichment <- tryCatch(
458
  )
459
 
460
  saveRDS(dds, file.path(outdir, "dds_fitted.rds"))
 
461
  summary <- list(
462
  input_kind = input_kind,
463
  normalization_method = if (is.null(txi)) "median_of_ratios" else "tximport_length_offsets",
 
54
  if (any(!vapply(covariates, safe_name, logical(1)))) stop("covariate name has unsupported characters")
55
  }
56
 
57
+ progress_path <- Sys.getenv("BIONEXUS_EXPRESSION_PROGRESS", "")
58
+ progress <- function(step, status) {
59
+ if (nzchar(progress_path)) cat(toJSON(list(step = as.integer(step), status = status), auto_unbox = TRUE), "\n", file = progress_path, append = TRUE, sep = "")
60
+ }
61
+ progress(1, "running")
62
  txi <- NULL
63
  if (input_kind == "salmon") {
64
  suppressPackageStartupMessages(library(tximport))
 
90
  }
91
  if (anyDuplicated(colnames(counts)) > 0) stop("count matrix contains duplicate sample names")
92
 
93
+ progress(1, "complete")
94
+ progress(2, "running")
95
  meta <- read.delim(metadata_path, check.names = FALSE, stringsAsFactors = FALSE)
96
  if (!("sample" %in% names(meta))) stop("metadata must contain a 'sample' column")
97
  if (!(condition_col %in% names(meta))) stop(paste0("metadata is missing condition column '", condition_col, "'"))
 
248
  if (length(design_warnings) == 0) design_audit$status <- "PASS"
249
  write(toJSON(design_audit, auto_unbox = TRUE, pretty = TRUE, null = "null", digits = 10), file.path(outdir, "design_audit.json"))
250
 
251
+ progress(2, "complete")
252
+ progress(3, "running")
253
  dds <- if (is.null(txi)) {
254
  DESeqDataSetFromMatrix(countData = counts, colData = meta, design = design_formula)
255
  } else {
256
  DESeqDataSetFromTximport(txi, colData = meta, design = design_formula)
257
  }
258
+ progress(3, "complete")
259
+ progress(4, "running")
260
  genes_input <- nrow(dds)
261
  keep <- rowSums(counts(dds) >= min_count) >= min_samples
262
  dds <- dds[keep, ]
263
  genes_kept <- nrow(dds)
264
  if (genes_kept < 2) stop("pre-filtering retained fewer than two genes")
265
 
266
+ progress(4, "complete")
267
+ progress(5, "running")
268
  dds <- estimateSizeFactors(dds)
269
  normalized <- counts(dds, normalized = TRUE)
270
  size_factors <- sizeFactors(dds)
 
288
  write.table(size_factor_df, file.path(outdir, "size_factors.tsv"), sep = "\t", quote = FALSE, row.names = FALSE)
289
  write.table(data.frame(gene = rownames(normalized), normalized, check.names = FALSE), file.path(outdir, "normalized_counts.tsv"), sep = "\t", quote = FALSE, row.names = FALSE)
290
 
291
+ progress(5, "complete")
292
+ progress(6, "running")
293
  # vst() is fast for ordinary whole-transcriptome matrices but requires enough rows
294
  # for its trend subsampling. varianceStabilizingTransformation() is the exact
295
  # DESeq2 transform and is robust for small deterministic validation subsets.
 
340
  quit(save = "no", status = 0)
341
  }
342
 
343
+ progress(6, "complete")
344
+ progress(7, "running")
345
  dds <- DESeq(dds)
346
  res <- results(dds, contrast = c(condition_col, test_level, reference_level), alpha = alpha)
347
  res_df <- as.data.frame(res)
 
375
  write.table(deg_df, file.path(outdir, "deseq2_significant.tsv"), sep = "\t", quote = FALSE, row.names = FALSE, na = "")
376
 
377
 
378
+ progress(7, "complete")
379
+ progress(8, "running")
380
  pca$.condition <- pca[[condition_col]]
381
  pca_plot <- ggplot(pca, aes(x = PC1, y = PC2, label = name, shape = .condition)) +
382
  geom_point(size = 3.2) +
 
418
  theme_minimal(base_size = 11) + theme(legend.position = "bottom")
419
  plot_svg_pdf_png("volcano", 8, 6, function() print(volcano_plot))
420
 
421
+ progress(8, "complete")
422
+ progress(9, "running")
423
  heatmap_written <- FALSE
424
  heatmap_basis <- "none"
425
  heatmap_gene_count <- 0L
 
466
  heatmap_written <- TRUE
467
  }
468
 
469
+ progress(9, if (heatmap_written) "complete" else "unavailable")
470
+ progress(10, "running")
471
  script_arg <- grep("^--file=", commandArgs(), value = TRUE)[[1]]
472
  source(file.path(dirname(sub("^--file=", "", script_arg)), "go_enrichment.R"))
473
  enrichment <- tryCatch(
 
481
  )
482
 
483
  saveRDS(dds, file.path(outdir, "dds_fitted.rds"))
484
+ progress(10, if (enrichment$status == "SUCCEEDED") "complete" else "unavailable")
485
  summary <- list(
486
  input_kind = input_kind,
487
  normalization_method = if (is.null(txi)) "median_of_ratios" else "tximport_length_offsets",
app/rnaseq/durable.py CHANGED
@@ -131,9 +131,31 @@ def status(job_id: str, user: str) -> dict:
131
  result["checkpoint_sha256"] = digest
132
  return {"job_id": job_id, "state": row["state"], "phase": row["phase"],
133
  "attempts": row["attempts"], "error": row["error"], "result": result,
 
134
  **_monitor(row)}
135
 
136
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
137
  def _monitor(row) -> dict:
138
  age = max(0.0, time.time() - row["heartbeat"]) if row["heartbeat"] else None
139
  return {"created_at": row["created"], "updated_at": row["updated"],
@@ -234,7 +256,8 @@ def run_job(job: dict) -> None:
234
  params=ExpressionParameters(**payload["params"]), source_label="reviewed-recovery",
235
  source_metadata={**payload["source"], "input_checksums": payload["checksums"], "qc_review": payload.get("qc_review")},
236
  stage=job["phase"], checkpoint_path=checkpoint, input_kind=payload["input_kind"],
237
- tx2gene_path=folder / "tx2gene.tsv" if payload["input_kind"] == "salmon" else None)
 
238
  if job["phase"] == "qc":
239
  result["checkpoint_sha256"] = sha256_file(checkpoint)
240
  result["checkpoint_file"] = checkpoint.name
 
131
  result["checkpoint_sha256"] = digest
132
  return {"job_id": job_id, "state": row["state"], "phase": row["phase"],
133
  "attempts": row["attempts"], "error": row["error"], "result": result,
134
+ "progress": read_progress(root() / job_id / f"progress_{row['attempts']}.jsonl"),
135
  **_monitor(row)}
136
 
137
 
138
+ def read_progress(path: Path) -> list[dict]:
139
+ """Only complete, bounded engine events from the current attempt are public."""
140
+ if not path.is_file():
141
+ return []
142
+ events = []
143
+ with path.open("rb") as stream:
144
+ body = stream.read(65536)
145
+ for line in body.splitlines(keepends=True):
146
+ if not line.endswith(b"\n"):
147
+ continue
148
+ try:
149
+ event = json.loads(line)
150
+ if (isinstance(event, dict) and type(event.get("step")) is int
151
+ and 1 <= event["step"] <= 10
152
+ and event.get("status") in {"running", "complete", "unavailable", "failed"}):
153
+ events.append({"step": event["step"], "status": event["status"]})
154
+ except (ValueError, UnicodeDecodeError):
155
+ continue
156
+ return events[-64:]
157
+
158
+
159
  def _monitor(row) -> dict:
160
  age = max(0.0, time.time() - row["heartbeat"]) if row["heartbeat"] else None
161
  return {"created_at": row["created"], "updated_at": row["updated"],
 
256
  params=ExpressionParameters(**payload["params"]), source_label="reviewed-recovery",
257
  source_metadata={**payload["source"], "input_checksums": payload["checksums"], "qc_review": payload.get("qc_review")},
258
  stage=job["phase"], checkpoint_path=checkpoint, input_kind=payload["input_kind"],
259
+ tx2gene_path=folder / "tx2gene.tsv" if payload["input_kind"] == "salmon" else None,
260
+ progress_path=folder / f"progress_{job['attempts'] + 1}.jsonl")
261
  if job["phase"] == "qc":
262
  result["checkpoint_sha256"] = sha256_file(checkpoint)
263
  result["checkpoint_file"] = checkpoint.name
app/rnaseq/expression.py CHANGED
@@ -236,7 +236,7 @@ def validate_artifacts(outdir: Path, summary: dict, stage: str) -> None:
236
  raise RnaSeqExpressionError("Required output artifacts missing: " + ", ".join(missing))
237
 
238
 
239
- def _run_r(counts_path: Path, metadata_path: Path, outdir: Path, params: ExpressionParameters, timeout_seconds: int = 600, *, stage: str = "full", checkpoint_path: Path | None = None, input_kind: str = "raw_counts", tx2gene_path: Path | None = None) -> dict[str, Any]:
240
  rscript = shutil.which(os.environ.get("BIONEXUS_RSCRIPT", "Rscript"))
241
  if not rscript:
242
  raise RnaSeqExpressionError("Rscript is not installed in the backend runtime; DESeq2 execution is unavailable.")
@@ -261,7 +261,8 @@ def _run_r(counts_path: Path, metadata_path: Path, outdir: Path, params: Express
261
  ]
262
  try:
263
  environment = dict(os.environ, BIONEXUS_EXPRESSION_STAGE=stage, BIONEXUS_EXPRESSION_INPUT=input_kind,
264
- BIONEXUS_EXPRESSION_CHECKPOINT=str(checkpoint_path or ""), BIONEXUS_TX2GENE=str(tx2gene_path or ""))
 
265
  completed = subprocess.run(argv, env=environment, capture_output=True, text=True, timeout=timeout_seconds, check=False)
266
  except subprocess.TimeoutExpired as exc:
267
  raise RnaSeqExpressionError(f"DESeq2 execution exceeded {timeout_seconds} seconds.") from exc
@@ -294,6 +295,7 @@ def execute_expression_analysis(
294
  checkpoint_path: Path | None = None,
295
  input_kind: str = "raw_counts",
296
  tx2gene_path: Path | None = None,
 
297
  ) -> dict[str, Any]:
298
  params.validate()
299
  if counts_path.stat().st_size > MAX_COUNTS_BYTES:
@@ -314,7 +316,8 @@ def execute_expression_analysis(
314
  outdir = Path(tmp) / "results"
315
  outdir.mkdir(parents=True, exist_ok=True)
316
  summary = _run_r(counts_path, metadata_path, outdir, params, stage=stage, checkpoint_path=checkpoint_path,
317
- input_kind=input_kind, tx2gene_path=tx2gene_path)
 
318
  validate_artifacts(outdir, summary, stage)
319
 
320
  provenance = {
 
236
  raise RnaSeqExpressionError("Required output artifacts missing: " + ", ".join(missing))
237
 
238
 
239
+ def _run_r(counts_path: Path, metadata_path: Path, outdir: Path, params: ExpressionParameters, timeout_seconds: int = 600, *, stage: str = "full", checkpoint_path: Path | None = None, input_kind: str = "raw_counts", tx2gene_path: Path | None = None, progress_path: Path | None = None) -> dict[str, Any]:
240
  rscript = shutil.which(os.environ.get("BIONEXUS_RSCRIPT", "Rscript"))
241
  if not rscript:
242
  raise RnaSeqExpressionError("Rscript is not installed in the backend runtime; DESeq2 execution is unavailable.")
 
261
  ]
262
  try:
263
  environment = dict(os.environ, BIONEXUS_EXPRESSION_STAGE=stage, BIONEXUS_EXPRESSION_INPUT=input_kind,
264
+ BIONEXUS_EXPRESSION_CHECKPOINT=str(checkpoint_path or ""), BIONEXUS_TX2GENE=str(tx2gene_path or ""),
265
+ BIONEXUS_EXPRESSION_PROGRESS=str(progress_path or ""))
266
  completed = subprocess.run(argv, env=environment, capture_output=True, text=True, timeout=timeout_seconds, check=False)
267
  except subprocess.TimeoutExpired as exc:
268
  raise RnaSeqExpressionError(f"DESeq2 execution exceeded {timeout_seconds} seconds.") from exc
 
295
  checkpoint_path: Path | None = None,
296
  input_kind: str = "raw_counts",
297
  tx2gene_path: Path | None = None,
298
+ progress_path: Path | None = None,
299
  ) -> dict[str, Any]:
300
  params.validate()
301
  if counts_path.stat().st_size > MAX_COUNTS_BYTES:
 
316
  outdir = Path(tmp) / "results"
317
  outdir.mkdir(parents=True, exist_ok=True)
318
  summary = _run_r(counts_path, metadata_path, outdir, params, stage=stage, checkpoint_path=checkpoint_path,
319
+ input_kind=input_kind, tx2gene_path=tx2gene_path,
320
+ **({"progress_path": progress_path} if progress_path else {}))
321
  validate_artifacts(outdir, summary, stage)
322
 
323
  provenance = {
app/routers/geo_search.py CHANGED
@@ -59,6 +59,7 @@ class AutomaticSelection(BaseModel):
59
  group_key: str | None = None
60
  source_sha256: str | None = None
61
  source_record_sha256: str | None = None
 
62
 
63
 
64
  async def _prepare_automatic(selection: AutomaticSelection):
@@ -96,6 +97,10 @@ async def prepare_geo_analysis(selection: AutomaticSelection):
96
  @router.post('/run')
97
  async def run_geo_analysis(selection: AutomaticSelection, user_id: str = Depends(require_user_id)):
98
  try:
 
 
 
 
99
  prepared, matrix = await _prepare_automatic(selection)
100
  if prepared['status'] != 'ready':
101
  raise GeoCountsError('Load the study and choose the indicated comparison before running.')
 
59
  group_key: str | None = None
60
  source_sha256: str | None = None
61
  source_record_sha256: str | None = None
62
+ require_durable: bool = False
63
 
64
 
65
  async def _prepare_automatic(selection: AutomaticSelection):
 
97
  @router.post('/run')
98
  async def run_geo_analysis(selection: AutomaticSelection, user_id: str = Depends(require_user_id)):
99
  try:
100
+ if selection.require_durable:
101
+ queue_status = await run_in_threadpool(durable.readiness)
102
+ if not queue_status['available']:
103
+ raise RnaSeqExpressionError('Live queued analysis unavailable: ' + ', '.join(queue_status['missing']))
104
  prepared, matrix = await _prepare_automatic(selection)
105
  if prepared['status'] != 'ready':
106
  raise GeoCountsError('Load the study and choose the indicated comparison before running.')
app/routers/jobs.py CHANGED
@@ -72,7 +72,7 @@ async def get_job(job_id: str, user_id: str = Depends(require_user_id)):
72
 
73
  # Hydrate from Storage if result was offloaded. Ownership has already been
74
  # established above; an unowned job never reaches artifact retrieval.
75
- if job.get("storage_url") and not job.get("context_json", {}).get("blast"):
76
  from app.services.artifact_storage import download_json
77
  results = download_json(job["storage_url"])
78
  if results:
 
72
 
73
  # Hydrate from Storage if result was offloaded. Ownership has already been
74
  # established above; an unowned job never reaches artifact retrieval.
75
+ if job.get("storage_url") and job.get("status") in {"complete", "failed"}:
76
  from app.services.artifact_storage import download_json
77
  results = download_json(job["storage_url"])
78
  if results:
app/routers/pipeline_v2.py CHANGED
@@ -82,6 +82,25 @@ def _is_real_uuid(job_id: str) -> bool:
82
  return False
83
 
84
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
85
  def _persist_v2_final(job_id: str, status: str, context: dict, error: str | None = None):
86
  """Update the persisted wizard job row with the final status + full context."""
87
  if not _is_real_uuid(job_id):
@@ -515,6 +534,8 @@ async def _execute(job_id: str, sequence: str, steps: list[str], status_callback
515
  _failed_step = "blast"
516
  _failed_error = result.get("error", "BLAST completed without usable hits")
517
  context["blast"] = result
 
 
518
 
519
  if denovo_mode:
520
  context["query"]["confidence"] = "de_novo"
 
82
  return False
83
 
84
 
85
+ def _persist_blast_preview(job_id: str, context: dict) -> None:
86
+ """Publish completed search evidence without declaring downstream work done."""
87
+ if not _is_real_uuid(job_id):
88
+ return
89
+ try:
90
+ client = get_supabase()
91
+ rows = client.table("jobs").select("user_id,context_json").eq("id", job_id).execute().data or []
92
+ if not rows:
93
+ return
94
+ stored = rows[0].get("context_json") or {}
95
+ if isinstance(stored, str):
96
+ import json
97
+ stored = json.loads(stored)
98
+ preview = {**(stored if isinstance(stored, dict) else {}), **context}
99
+ client.table("jobs").update({"context_json": preview}).eq("id", job_id).eq("user_id", rows[0]["user_id"]).execute()
100
+ except Exception:
101
+ logger.warning("Could not persist early BLAST results for %s", job_id)
102
+
103
+
104
  def _persist_v2_final(job_id: str, status: str, context: dict, error: str | None = None):
105
  """Update the persisted wizard job row with the final status + full context."""
106
  if not _is_real_uuid(job_id):
 
534
  _failed_step = "blast"
535
  _failed_error = result.get("error", "BLAST completed without usable hits")
536
  context["blast"] = result
537
+ if not provider_failed:
538
+ await asyncio.to_thread(_persist_blast_preview, job_id, context)
539
 
540
  if denovo_mode:
541
  context["query"]["confidence"] = "de_novo"