metabolomics: run targeted panels in raw mode (platforms, internal standards, LOD, tank units)

.github/workflows/ci.yml

@@ -141,6 +141,14 @@ jobs:
             any::knitr
             any::rmarkdown
 
+      # Not through setup-r-dependencies: pak cannot solve bioc::limma there
+      # (limma suggests vsn, which imports limma). Hard dependencies only.
+      - name: Install limma
+        run: |
+          Rscript -e 'install.packages("BiocManager", repos = "https://cloud.r-project.org")' \
+                  -e 'BiocManager::install("limma", update = FALSE, ask = FALSE, dependencies = c("Depends", "Imports", "LinkingTo"))' \
+                  -e 'stopifnot(requireNamespace("limma", quietly = TRUE)); cat("limma", as.character(packageVersion("limma")), "\n")'
+
       - uses: quarto-dev/quarto-actions/setup@v2
 
       - name: Install current Nextflow
@@ -220,6 +228,24 @@ jobs:
           test "$(ls /tmp/aree-smoke/rnaseq-raw-se/rnaseq/trimmed/*_trimmed_R1.fastq.gz | wc -l)" -eq 4
           ! ls /tmp/aree-smoke/rnaseq-raw-se/rnaseq/trimmed/*_trimmed_R2.fastq.gz 2>/dev/null
 
+      # A real (not stub) run: normalization, the detection filter, tank
+      # units and limma all execute on the simulated targeted-panel fixture.
+      - name: Run the metabolomics raw DAG on a targeted-panel fixture
+        run: |
+          nextflow run workflows/metabolomics -profile local \
+            --mode raw_reanalysis \
+            --metabolomics.raw_feature_table tests/fixtures/metabolomics_raw/features.tsv \
+            --metabolomics.sample_sheet tests/fixtures/metabolomics_raw/samplesheet.tsv \
+            --metabolomics.study_id CI_METAB \
+            --metabolomics.comparison_id treatment_vs_control \
+            --metabolomics.normalization internal_standard \
+            --metabolomics.detection_limit 0.001 \
+            --metabolomics.replicate_unit tank \
+            --metabolomics.acquisition targeted \
+            --outdir /tmp/aree-smoke/metabolomics-raw \
+            -work-dir /tmp/aree-smoke/work-metabolomics-raw
+          python scripts/check_metabolomics_raw_smoke.py /tmp/aree-smoke/metabolomics-raw
+
       # Two stub runs differing only in comparison params, paired and single
       # end: the second must take FASTQC, fastp, Salmon and viral quant from
       # the -resume cache. Skipped in the Python job, which has no Nextflow.
@@ -258,6 +284,9 @@ jobs:
           python scripts/check_manifest_versions.py \
             /tmp/aree-smoke/rnaseq-raw/rnaseq/manifests/CI_RNASEQ_treatment_vs_control_manifest.json \
             /tmp/aree-smoke/work-rnaseq-raw
+          python scripts/check_manifest_versions.py \
+            /tmp/aree-smoke/metabolomics-raw/metabolomics/manifest/CI_METAB_treatment_vs_control_manifest.json \
+            /tmp/aree-smoke/work-metabolomics-raw
 
   schema-validation:
     name: JSON schema sanity

docs/implementation_status.md

@@ -44,7 +44,7 @@ demo pipeline runs from the committed header-only registry).
 | RNA-seq Nextflow workflow | `workflows/rnaseq/`, `modules/rnaseq/` | **Executed at full depth against real public FASTQ** on eight comparisons (three from PRJNA1329250, one from PRJNA593309, three from PRJNA826964 that failed QC, one from PRJNA735889): FASTQC → fastp → Salmon quant → MultiQC → DESeq2 → standardize → viral QC → manifest → report, 23,189–30,625 genes per comparison, every one with a real `lfcSE`. An arm can pool several condition levels, and `--replicate_unit` fits DESeq2 on pseudo-bulk units (e.g. tanks) where treatment was applied per unit; the manifest records both. The two 2026-09-30 CALLA2026 re-runs are the first to finish `OK`, with the report task in local scratch and the work directory off exFAT; every earlier full-depth run exited `ERR` at `RENDER_REPORT`. Processed mode runs end to end; CI also stub-runs the raw DAG | runs are native (`-profile local`, Homebrew tools), **not** inside the declared containers; manifests committed before 2026-10-11 list the declared container versions in `software_versions` (e.g. Salmon 1.10.3, DESeq2 1.42.0), not those that ran (Salmon 2.6.0, DESeq2 1.50.2, tximport 1.38.2); newer runs build it from each process's `versions.yml`, and CI checks the manifest against the work dir (`scripts/check_manifest_versions.py`); `qc_metrics` is still null; a non-stub run on the tiny CI fixtures fails in `SALMON_QUANT` under Salmon 2.6; the CI stub run validates wiring, not tools or biology |
 | Methylation Nextflow workflow | `workflows/methylation/`, `modules/methylation/` | parses and runs `processed_results_harmonization` end to end; real Bismark/methylKit command lines | raw path never run against bisulfite data; R scripts unexecuted |
 | Proteomics Nextflow workflow | `workflows/proteomics/`, `modules/proteomics/` | parses and runs `processed_results_harmonization` end to end; limma diff-abundance written | raw path not run; no raw abundance-matrix fixture |
-| Metabolomics Nextflow workflow | `workflows/metabolomics/`, `modules/metabolomics/` | parses and runs `processed_results_harmonization` end to end | raw path not run against real feature-table data |
+| Metabolomics Nextflow workflow | `workflows/metabolomics/`, `modules/metabolomics/` | runs `processed_results_harmonization` end to end; CI runs the raw path for real (not stubbed) on a simulated targeted-panel fixture: per-platform internal-standard normalization, detection limit, `--replicate_unit`, limma with `lfcSE`, missingness report, manifest and report | native only (`-profile local`), no container pulled; the untargeted TIC path has no real-data run |
 | Container images | `containers/README.md` | declared image tags per step | **no image has ever been pulled or built.** Every real execution used Homebrew-installed tools via `-profile local`, so the container specifications remain entirely unverified |
 
 `processed_results_harmonization` is now **verified** for all four workflows —
@@ -88,9 +88,10 @@ worth less than one that says what it got wrong.
   `resolve_identifier` so that legacy symbols from older publications resolve.
 - Additional species beyond *C. gigas* (schema already carries `species`; this
   is a data + reference addition, see [adding_a_species.md](adding_a_species.md)).
-- Real-data execution of the methylation, proteomics, and metabolomics raw
-  paths. RNA-seq has tiny paired-FASTQ fixtures for a stubbed CI DAG test;
-  equivalent fixtures for the other assays remain future work.
+- Real-data execution of the methylation and proteomics raw paths. RNA-seq
+  has tiny paired-FASTQ fixtures for a stubbed CI DAG test and metabolomics a
+  targeted-panel fixture CI runs for real; equivalent fixtures for methylation
+  and proteomics remain future work.
 - Recording observed tool versions and parsed Salmon/MultiQC QC metrics in the
   RNA-seq manifest, instead of declared versions and null `qc_metrics`.
 - A covariance-aware within-study model, so a study can contribute more than

modules/metabolomics/differential_abundance.nf

@@ -4,12 +4,18 @@
 // small-n omics intensity matrices, widely used for metabolomics feature
 // tables as well as microarray/proteomics data.
 //
-// STATUS: structurally complete DSL2 process. The R script below is real,
-// syntactically valid limma usage (lmFit -> eBayes -> topTable) that would
-// run correctly given a real normalized intensity matrix and sample sheet.
-// NOT executed in this build — no compute / no real normalized-matrix input
-// exists here (see workflows/metabolomics/README.md).
-
+// With --metabolomics.replicate_unit (a sample-sheet column, e.g. tank), log2
+// values are averaged within each unit before the fit, so the units are the
+// replicates; a unit may not sit in both arms. Missing values (below
+// detection) are left out of each unit's mean. A feature is tested only if it
+// is detected in at least --metabolomics.min_detected_per_arm replicates in
+// each arm; every feature, tested or not, is listed in the missingness report.
+//
+// lfcSE is limma's posterior standard error of log2FC
+// (sqrt(s2.post) * stdev.unscaled), so the evidence record carries one.
+//
+// STATUS: executed on the real targeted panel of DUPERRET2025_POMS_STARVE
+// (2026-10-16) and by the raw-mode CI smoke test.
 process DIFFERENTIAL_ABUNDANCE {
     tag "${study_id}:${comparison_id}"
     label 'process_medium'
@@ -21,83 +27,145 @@ process DIFFERENTIAL_ABUNDANCE {
     val comparison_id
     path normalized_table
     path sample_sheet
+    val replicate_unit
+    val min_detected_per_arm
+    val acquisition
 
     output:
     tuple val(study_id), val(comparison_id), path("${study_id}_${comparison_id}_diffabundance_raw.tsv"), emit: results
+    path "${study_id}_${comparison_id}_missingness.tsv", emit: missingness
+    path "${study_id}_${comparison_id}_diffabundance_qc.json", emit: da_qc
     path "versions.yml", emit: versions
 
     script:
-    // Real limma moderated-t-test skeleton. Output columns are written to
-    // already match the STANDARDIZE_OUTPUT target schema (feature_id,
-    // log2FC, pvalue, padj) plus metadata columns carried through so that
-    // step is largely a pass-through/validation in raw mode.
     """
     cat <<-'EOF' > run_limma.R
     #!/usr/bin/env Rscript
     # AREE metabolomics differential abundance (limma moderated t-test)
-    #
-    # NOTE: structurally complete, unexecuted in this build. Written to run
-    # correctly against a real log2-normalized feature x sample intensity
-    # matrix and a real two-group sample sheet.
     suppressPackageStartupMessages({
         library(limma)
+        library(jsonlite)
     })
 
     args <- list(
         normalized_table = "${normalized_table}",
         sample_sheet     = "${sample_sheet}",
         study_id         = "${study_id}",
         comparison_id    = "${comparison_id}",
-        out_tsv          = "${study_id}_${comparison_id}_diffabundance_raw.tsv"
+        replicate_unit   = "${replicate_unit ?: ''}",
+        min_detected     = as.integer("${min_detected_per_arm}"),
+        acquisition      = "${acquisition}",
+        out_tsv          = "${study_id}_${comparison_id}_diffabundance_raw.tsv",
+        out_missing      = "${study_id}_${comparison_id}_missingness.tsv",
+        out_qc           = "${study_id}_${comparison_id}_diffabundance_qc.json"
     )
+    if (!(args\$acquisition %in% c("untargeted", "targeted"))) {
+        stop("acquisition must be 'untargeted' or 'targeted', got '", args\$acquisition, "'")
+    }
 
     features <- read.delim(args\$normalized_table, stringsAsFactors = FALSE, check.names = FALSE)
-    samples  <- read.delim(args\$sample_sheet, stringsAsFactors = FALSE)
+    samples  <- read.delim(args\$sample_sheet, stringsAsFactors = FALSE, check.names = FALSE)
     samples\$condition <- factor(samples\$condition, levels = c("control", "treatment"))
+    if (any(is.na(samples\$condition))) {
+        stop("sample sheet condition must be 'control' or 'treatment' for every sample")
+    }
 
     meta_cols <- intersect(
-        c("feature_id", "putative_metabolite_name", "annotation_confidence_level"),
+        c("feature_id", "putative_metabolite_name", "annotation_confidence_level", "platform"),
         colnames(features)
     )
     sample_ids <- intersect(samples\$sample_id, colnames(features))
     stopifnot(length(sample_ids) >= 2)
-
     samples <- samples[match(sample_ids, samples\$sample_id), ]
-
     expr <- as.matrix(features[, sample_ids, drop = FALSE])
     rownames(expr) <- features\$feature_id
 
-    design <- model.matrix(~condition, data = samples)
-
-    fit <- lmFit(expr, design)
-    fit <- eBayes(fit)
-
-    tt <- topTable(fit, coef = "conditiontreatment", number = Inf, sort.by = "none")
+    unit_col <- trimws(args\$replicate_unit)
+    if (nzchar(unit_col)) {
+        if (!(unit_col %in% colnames(samples))) {
+            stop("replicate_unit column '", unit_col, "' is not in the sample sheet")
+        }
+        units <- as.character(samples[[unit_col]])
+        if (any(is.na(units) | units == "")) {
+            stop("replicate_unit column '", unit_col, "' is empty for some samples")
+        }
+        arms_per_unit <- tapply(as.character(samples\$condition), units, function(x) length(unique(x)))
+        if (any(arms_per_unit > 1)) {
+            stop("replicate_unit(s) spanning both arms: ",
+                 paste(names(arms_per_unit)[arms_per_unit > 1], collapse = ", "))
+        }
+        unit_ids <- unique(units)
+        expr <- sapply(unit_ids, function(u) {
+            m <- rowMeans(expr[, units == u, drop = FALSE], na.rm = TRUE)
+            m[is.nan(m)] <- NA
+            m
+        })
+        expr <- matrix(expr, nrow = nrow(features), dimnames = list(features\$feature_id, unit_ids))
+        condition <- factor(as.character(samples\$condition[match(unit_ids, units)]),
+                            levels = c("control", "treatment"))
+    } else {
+        condition <- samples\$condition
+    }
 
+    detected <- !is.na(expr)
+    n_control   <- rowSums(detected[, condition == "control", drop = FALSE])
+    n_treatment <- rowSums(detected[, condition == "treatment", drop = FALSE])
+    tested <- n_control >= args\$min_detected & n_treatment >= args\$min_detected
+    reason <- ifelse(tested, "tested",
+              ifelse(n_control == 0 & n_treatment == 0, "not_detected",
+              ifelse(n_control < args\$min_detected & n_treatment < args\$min_detected,
+                     "detected_too_rarely_in_both_arms",
+              ifelse(n_control < args\$min_detected, "detected_too_rarely_in_control",
+                     "detected_too_rarely_in_treatment"))))
+    missingness <- data.frame(features[, meta_cols, drop = FALSE],
+                              n_replicates_control = sum(condition == "control"),
+                              n_detected_control = n_control,
+                              n_replicates_treatment = sum(condition == "treatment"),
+                              n_detected_treatment = n_treatment,
+                              status = reason, check.names = FALSE, stringsAsFactors = FALSE)
+    write.table(missingness, args\$out_missing, sep = "\\t", quote = FALSE, row.names = FALSE, na = "NA")
+    if (sum(tested) < 2) stop("fewer than 2 features pass min_detected_per_arm = ", args\$min_detected)
+
+    design <- model.matrix(~condition)
+    fit <- eBayes(lmFit(expr[tested, , drop = FALSE], design))
+    coef <- "conditiontreatment"
+    tt <- topTable(fit, coef = coef, number = Inf, sort.by = "none")
     out <- data.frame(
-        feature_id = rownames(tt),
-        log2FC     = tt\$logFC,
-        pvalue     = tt\$P.Value,
-        padj       = tt\$adj.P.Val,
+        feature_id      = rownames(tt),
+        log2FC          = tt\$logFC,
+        lfcSE           = sqrt(fit\$s2.post) * fit\$stdev.unscaled[, coef],
+        pvalue          = tt\$P.Value,
+        padj            = tt\$adj.P.Val,
+        analysis_method = paste0(args\$acquisition, "_MS_limma_moderated_t"),
         stringsAsFactors = FALSE
     )
-
     # Carry through annotation metadata already attached upstream so
     # downstream steps do not need to re-join it.
-    if ("putative_metabolite_name" %in% meta_cols) {
-        out\$putative_metabolite_name <- features\$putative_metabolite_name[match(out\$feature_id, features\$feature_id)]
-    }
-    if ("annotation_confidence_level" %in% meta_cols) {
-        out\$annotation_confidence_level <- features\$annotation_confidence_level[match(out\$feature_id, features\$feature_id)]
+    for (col in setdiff(meta_cols, "feature_id")) {
+        out[[col]] <- features[[col]][match(out\$feature_id, features\$feature_id)]
     }
-
     out <- out[order(out\$padj, na.last = TRUE), ]
-
     write.table(out, args\$out_tsv, sep = "\\t", quote = FALSE, row.names = FALSE, na = "NA")
 
+    qc <- list(
+        study_id = args\$study_id,
+        comparison_id = args\$comparison_id,
+        replicate_unit = if (nzchar(unit_col)) unit_col else NULL,
+        n_replicates_control = sum(condition == "control"),
+        n_replicates_treatment = sum(condition == "treatment"),
+        min_detected_per_arm = args\$min_detected,
+        n_features_in = nrow(features),
+        n_features_tested = sum(tested),
+        n_features_not_tested = as.list(table(reason[!tested])),
+        n_padj_below_0.05 = sum(out\$padj < 0.05, na.rm = TRUE),
+        limma_df_prior = unname(fit\$df.prior[1]),
+        analysis_method = out\$analysis_method[1]
+    )
+    writeLines(toJSON(qc, auto_unbox = TRUE, pretty = TRUE, null = "null"), args\$out_qc)
     cat(sprintf(
-        "AREE limma differential abundance complete for %s:%s — %d features tested\\n",
-        args\$study_id, args\$comparison_id, nrow(out)
+        "AREE limma differential abundance complete for %s:%s: %d of %d features tested on %d vs %d replicates\\n",
+        args\$study_id, args\$comparison_id, sum(tested), nrow(features),
+        sum(condition == "control"), sum(condition == "treatment")
     ))
     EOF
 
@@ -107,12 +175,15 @@ process DIFFERENTIAL_ABUNDANCE {
     ${task.process}:
         r-base: \$(Rscript -e 'cat(as.character(getRversion()))')
         bioconductor-limma: \$(Rscript -e 'cat(as.character(packageVersion("limma")))')
+        r-jsonlite: \$(Rscript -e 'cat(as.character(packageVersion("jsonlite")))')
     END_VERSIONS
     """
 
     stub:
     """
     touch ${study_id}_${comparison_id}_diffabundance_raw.tsv
+    touch ${study_id}_${comparison_id}_missingness.tsv
+    touch ${study_id}_${comparison_id}_diffabundance_qc.json
     touch versions.yml
     """
 }

modules/metabolomics/emit_manifest.nf

@@ -149,6 +149,9 @@ process EMIT_MANIFEST {
                              "separately by src/harmonize/metabolomics.py",
         "parameters": {
             "mode": MODE,
+            # Every params.metabolomics value this run used (normalization,
+            # detection_limit, replicate_unit, ...), not just the mode.
+            "metabolomics": json.loads(r'''${groovy.json.JsonOutput.toJson(params.metabolomics)}'''),
         },
         "input_checksums_sha256": checksums,
         # Tools as reported by each process that ran, keyed by process name.