diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 82ca9c2..481082e 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/docs/implementation_status.md b/docs/implementation_status.md index b292004..d9df2d3 100644 --- a/docs/implementation_status.md +++ b/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 diff --git a/modules/metabolomics/differential_abundance.nf b/modules/metabolomics/differential_abundance.nf index f606347..71f0aec 100644 --- a/modules/metabolomics/differential_abundance.nf +++ b/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,26 +27,24 @@ 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( @@ -48,56 +52,120 @@ process DIFFERENTIAL_ABUNDANCE { 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 """ } diff --git a/modules/metabolomics/emit_manifest.nf b/modules/metabolomics/emit_manifest.nf index 720ce05..13c48ae 100644 --- a/modules/metabolomics/emit_manifest.nf +++ b/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. diff --git a/modules/metabolomics/normalize_qc.nf b/modules/metabolomics/normalize_qc.nf index 9456949..d89502c 100644 --- a/modules/metabolomics/normalize_qc.nf +++ b/modules/metabolomics/normalize_qc.nf @@ -1,12 +1,26 @@ // Normalize feature intensities and compute QC metrics (raw_reanalysis mode // only). // -// Normalization: total-ion-current (TIC) normalization per sample (each -// sample's intensities are scaled so the sample's total intensity equals the -// across-sample median total intensity), followed by log2(x + 1) transform. -// This is a standard, simple untargeted-metabolomics normalization choice — -// not the only defensible one (median-fold-change / PQN are common -// alternatives) — documented here so a curator can swap it out per study. +// Normalization (--metabolomics.normalization): +// tic (default) each sample's intensities are scaled so its total equals +// the across-sample median total, then log2(x + 1). A +// standard, simple untargeted-metabolomics choice. +// median log2, then each sample is centred on its median +// detected feature (plus the grand median, to keep the +// scale). Suits targeted panels with no usable standard. +// internal_standard log2, then each sample's mean log2 internal-standard +// area is subtracted (plus its across-sample mean). The +// usual choice for targeted MRM panels that spike +// standards into every sample. +// Every method runs separately within each value of the optional `platform` +// column, so panels from different instruments are never scaled against each +// other. Internal-standard rows (`is_internal_standard` true) are dropped from +// the output: they are not metabolites to test. +// +// Detection limit (--metabolomics.detection_limit): values at or below it are +// below detection and become missing (NaN), not a tiny intensity. Without it, +// 'tic' keeps its original behaviour; 'median' and 'internal_standard' treat +// non-positive values as missing. // // QC: if the sample sheet marks any samples as QC pooled samples // (`qc_pool` column == true/1), computes the coefficient of variation (CV) @@ -16,7 +30,8 @@ // present, this metric is explicitly reported as not_applicable rather than // silently omitted. // -// STATUS: real, runnable pandas logic. NOT executed in this build. +// STATUS: 'tic' runs in no real study yet; 'internal_standard' with a +// detection limit is what DUPERRET2025_POMS_STARVE's targeted panel uses. process NORMALIZE_QC { tag "${study_id}:${comparison_id}" @@ -30,6 +45,8 @@ process NORMALIZE_QC { path annotated_table path sample_sheet val qc_cv_threshold + val normalization + val detection_limit output: tuple val(study_id), val(comparison_id), path("${study_id}_${comparison_id}_features_normalized.tsv"), emit: normalized_table @@ -41,10 +58,11 @@ process NORMALIZE_QC { cat <<-'EOF' > normalize_qc.py #!/usr/bin/env python3 ${'"""'} - TIC normalization + log2 transform of metabolite feature intensities, - plus pooled-QC-sample CV reporting. + Per-platform normalization + log2 transform of metabolite feature + intensities, plus pooled-QC-sample CV reporting. ${'"""'} import json + import sys import numpy as np import pandas as pd @@ -54,29 +72,72 @@ process NORMALIZE_QC { IN_PATH = "${annotated_table}" SAMPLE_SHEET_PATH = "${sample_sheet}" CV_THRESHOLD = float("${qc_cv_threshold}") + METHOD = "${normalization}" + DETECTION_LIMIT = None if "${detection_limit}" in ("", "null") else float("${detection_limit}") OUT_TABLE = "${study_id}_${comparison_id}_features_normalized.tsv" OUT_QC = "${study_id}_${comparison_id}_normalization_qc.json" - df = pd.read_csv(IN_PATH, sep="\\t") - samples = pd.read_csv(SAMPLE_SHEET_PATH, sep="\\t" if "\\t" in open(SAMPLE_SHEET_PATH).readline() else ",") + if METHOD not in ("tic", "median", "internal_standard"): + sys.exit(f"ERROR: normalization must be 'tic', 'median' or 'internal_standard', got {METHOD!r}") + + df = pd.read_csv(IN_PATH, sep="\t") + samples = pd.read_csv(SAMPLE_SHEET_PATH, sep="\t" if "\t" in open(SAMPLE_SHEET_PATH).readline() else ",") - meta_cols = [c for c in ["feature_id", "putative_metabolite_name", "annotation_confidence_level"] if c in df.columns] + meta_cols = [c for c in ["feature_id", "putative_metabolite_name", "annotation_confidence_level", "platform"] + if c in df.columns] sample_ids = samples["sample_id"].astype(str).tolist() intensity_cols = [c for c in sample_ids if c in df.columns] intensities = df[intensity_cols].astype(float) - - # Total-ion-current normalization: scale each sample column so its sum - # equals the median per-sample sum across all samples. - sample_totals = intensities.sum(axis=0) - median_total = sample_totals.median() - scaling_factors = median_total / sample_totals.replace(0, np.nan) - normalized = intensities.mul(scaling_factors, axis=1) - - log2_normalized = np.log2(normalized.clip(lower=0) + 1) - - out = pd.concat([df[meta_cols], log2_normalized], axis=1) - out.to_csv(OUT_TABLE, sep="\\t", index=False) + is_std = (df["is_internal_standard"].astype(str).str.lower().isin(["true", "1", "yes"]) + if "is_internal_standard" in df.columns else pd.Series(False, index=df.index)) + platform = df["platform"].astype(str) if "platform" in df.columns else pd.Series("all", index=df.index) + + n_below_detection = 0 + if DETECTION_LIMIT is not None: + below = intensities <= DETECTION_LIMIT + n_below_detection = int(below[~is_std].sum().sum()) + intensities = intensities.mask(below) + + normalized_log2 = pd.DataFrame(index=df.index, columns=intensity_cols, dtype=float) + normalized = pd.DataFrame(index=df.index, columns=intensity_cols, dtype=float) # linear, for the pooled-QC CV + scaling = {} + for plat in sorted(platform.unique()): + rows = platform == plat + feats = rows & ~is_std + block = intensities.loc[rows] + if METHOD == "tic": + totals = intensities.loc[feats].sum(axis=0) + factors = totals.median() / totals.replace(0, np.nan) + scaled = block.mul(factors, axis=1) + normalized.loc[rows] = scaled + if DETECTION_LIMIT is None: + normalized_log2.loc[rows] = np.log2(scaled.clip(lower=0) + 1) + else: + normalized_log2.loc[rows] = np.log2(scaled) + scaling[plat] = {k: (None if pd.isna(v) else round(float(v), 6)) for k, v in factors.items()} + continue + logged = np.log2(block.where(block > 0)) + if METHOD == "median": + ref = logged.loc[feats[rows]].median(axis=0) + else: + std_rows = is_std[rows] + if not std_rows.any(): + sys.exit(f"ERROR: normalization 'internal_standard' but platform {plat!r} has no " + "row with is_internal_standard true") + ref = logged.loc[std_rows].mean(axis=0) + if ref.isna().any(): + sys.exit(f"ERROR: platform {plat!r}: no detected internal standard in sample(s) " + f"{list(ref[ref.isna()].index)}") + offset = ref - ref.mean() + normalized_log2.loc[rows] = logged.sub(offset, axis=1) + normalized.loc[rows] = np.exp2(normalized_log2.loc[rows]) + scaling[plat] = {k: round(float(v), 6) for k, v in offset.items()} + + keep = ~is_std + log2_normalized = normalized_log2.loc[keep] + out = pd.concat([df.loc[keep, meta_cols], log2_normalized], axis=1) + out.to_csv(OUT_TABLE, sep="\t", index=False, na_rep="NA") qc_pool_col = samples.get("qc_pool") if qc_pool_col is not None: @@ -87,7 +148,7 @@ process NORMALIZE_QC { qc_pool_samples = [] if len(qc_pool_samples) >= 2: - pool_data = normalized[qc_pool_samples] # CV computed on linear scale, not log + pool_data = normalized.loc[keep, qc_pool_samples] # CV computed on linear scale, not log pool_mean = pool_data.mean(axis=1) pool_sd = pool_data.std(axis=1, ddof=1) cv = (pool_sd / pool_mean.replace(0, np.nan)).abs() @@ -111,17 +172,25 @@ process NORMALIZE_QC { qc = { "study_id": STUDY_ID, "comparison_id": COMPARISON_ID, - "normalization_method": "total_ion_current_then_log2", - "n_features": len(df), + "normalization_method": {"tic": "total_ion_current_then_log2", "median": "log2_then_sample_median", + "internal_standard": "log2_then_internal_standard"}[METHOD] + + ("_per_platform" if "platform" in df.columns else ""), + "detection_limit": DETECTION_LIMIT, + "n_values_below_detection": n_below_detection, + "n_features": int(keep.sum()), + "n_internal_standards_excluded": int(is_std.sum()), + "internal_standards": df.loc[is_std, "feature_id"].astype(str).tolist(), "n_samples_normalized": len(intensity_cols), - "sample_scaling_factors": {k: (None if pd.isna(v) else round(float(v), 4)) for k, v in scaling_factors.items()}, + # tic: linear scaling factor per sample; median/internal_standard: the + # log2 offset subtracted from each sample. Keyed by platform. + "sample_scaling": scaling, "qc_pool_metric": qc_pool_metric, } with open(OUT_QC, "w") as fh: json.dump(qc, fh, indent=2) print( - f"Normalized {len(intensity_cols)} samples x {len(df)} features for " + f"Normalized ({METHOD}) {len(intensity_cols)} samples x {int(keep.sum())} features for " f"{STUDY_ID}:{COMPARISON_ID}; QC pool metric: {qc_pool_metric}" ) EOF diff --git a/modules/metabolomics/standardize_output.nf b/modules/metabolomics/standardize_output.nf index eb5d8a2..1900a36 100644 --- a/modules/metabolomics/standardize_output.nf +++ b/modules/metabolomics/standardize_output.nf @@ -61,6 +61,9 @@ process STANDARDIZE_OUTPUT { "pvalue", "padj", ] + # Carried through when present: the standard error of log2FC (lets the + # evidence record be pooled) and the method that produced the row. + OPTIONAL = ["lfcSE", "analysis_method"] IN_PATH = "${raw_results}" OUT_PATH = "${study_id}_${comparison_id}_metabolite_features_standardized.tsv" @@ -76,6 +79,7 @@ process STANDARDIZE_OUTPUT { "log2FoldChange": "log2FC", "logFC": "log2FC", "p.value": "pvalue", "PValue": "pvalue", "p_value": "pvalue", "FDR": "padj", "adj.P.Val": "padj", "adjusted_p_value": "padj", + "standard_error": "lfcSE", "SE": "lfcSE", "se": "lfcSE", } df = df.rename(columns={k: v for k, v in rename_map.items() if k in df.columns}) @@ -87,7 +91,7 @@ process STANDARDIZE_OUTPUT { f"input. Present columns: {list(df.columns)}" ) - out = df[REQUIRED].copy() + out = df[REQUIRED + [c for c in OPTIONAL if c in df.columns]].copy() out["annotation_confidence_level"] = pd.to_numeric( out["annotation_confidence_level"], errors="coerce" @@ -105,7 +109,7 @@ process STANDARDIZE_OUTPUT { out.loc[bad_mask, "annotation_confidence_level"] = 4 out["annotation_confidence_level"] = out["annotation_confidence_level"].astype(int) - for col in ["log2FC", "pvalue", "padj"]: + for col in ["log2FC", "pvalue", "padj"] + (["lfcSE"] if "lfcSE" in out.columns else []): out[col] = pd.to_numeric(out[col], errors="coerce") out["metabolite_class"] = out["metabolite_class"].fillna("unknown") diff --git a/scripts/check_metabolomics_raw_smoke.py b/scripts/check_metabolomics_raw_smoke.py new file mode 100644 index 0000000..fe8cd89 --- /dev/null +++ b/scripts/check_metabolomics_raw_smoke.py @@ -0,0 +1,64 @@ +"""Check the metabolomics raw-mode CI run on tests/fixtures/metabolomics_raw. + +The fixture plants three effects (LC_up +2.0, LC_down -1.5, GC_up +1.0), one +metabolite detected only in treatment and one never detected. The run must +recover the effects in direction and rough size, give every tested row a +standard error, keep internal standards and undetected metabolites out of the +results, and list every non-standard feature in the missingness report. +""" +from __future__ import annotations + +import json +import sys +from pathlib import Path + +import pandas as pd + +PLANTED = {"LC_up": 2.0, "LC_down": -1.5, "GC_up": 1.0} +NOT_TESTED = {"LC_treatonly": "detected_too_rarely_in_control", "LC_lod": "not_detected"} +INTERNAL_STANDARDS = {"LC_IS", "GC_IS"} + + +def main(outdir: str) -> int: + root = Path(outdir) / "metabolomics" + stem = "CI_METAB_treatment_vs_control" + res = pd.read_csv(root / "standardized" / f"{stem}_metabolite_features_standardized.tsv", sep="\t") + miss = pd.read_csv(root / "differential_abundance" / f"{stem}_missingness.tsv", sep="\t") + manifest = json.loads((root / "manifest" / f"{stem}_manifest.json").read_text()) + errors = [] + + by_id = res.set_index("feature_id") + for fid, lfc in PLANTED.items(): + row = by_id.loc[fid] + if not (row.padj < 0.05 and abs(row.log2FC - lfc) < 0.5): + errors.append(f"{fid}: log2FC {row.log2FC:.2f}, padj {row.padj:.2g}; planted {lfc:+.1f}") + if res["lfcSE"].isna().any() or (res["lfcSE"] <= 0).any(): + errors.append("every tested row needs a positive lfcSE") + if set(res["analysis_method"]) != {"targeted_MS_limma_moderated_t"}: + errors.append(f"analysis_method: {sorted(set(res['analysis_method']))}") + leaked = (set(NOT_TESTED) | INTERNAL_STANDARDS) & set(res["feature_id"]) + if leaked: + errors.append(f"untestable features in the results: {sorted(leaked)}") + + status = miss.set_index("feature_id")["status"] + for fid, expected in NOT_TESTED.items(): + if status.get(fid) != expected: + errors.append(f"missingness {fid}: {status.get(fid)!r}, expected {expected!r}") + if INTERNAL_STANDARDS & set(miss["feature_id"]): + errors.append("internal standards belong in normalization QC, not the missingness report") + if set(miss["n_replicates_control"]) != {3}: + errors.append("replicate_unit tank should give 3 control units") + + params = manifest["parameters"]["metabolomics"] + for key, value in {"normalization": "internal_standard", "replicate_unit": "tank"}.items(): + if params.get(key) != value: + errors.append(f"manifest parameters.metabolomics.{key} = {params.get(key)!r}") + + for e in errors: + print(f"FAIL: {e}") + print(f"{len(res)} tested, {int((status != 'tested').sum())} not tested, {len(errors)} problem(s)") + return 1 if errors else 0 + + +if __name__ == "__main__": + sys.exit(main(sys.argv[1])) diff --git a/src/harmonize/metabolomics.py b/src/harmonize/metabolomics.py index 5f092ef..e195296 100644 --- a/src/harmonize/metabolomics.py +++ b/src/harmonize/metabolomics.py @@ -45,6 +45,10 @@ def harmonize_metabolomics( standardized = r.get("putative_metabolite_name") if confidence != "unresolved" else None effect_size = float(r["log2FC"]) if pd.notna(r.get("log2FC")) else None + standard_error = float(r["lfcSE"]) if pd.notna(r.get("lfcSE")) else None + # The workflow's raw mode writes the method per row; a processed table + # without the column falls back to the caller's default. + method = r.get("analysis_method") if pd.notna(r.get("analysis_method")) else analysis_method direction = molecular_direction_from_effect(effect_size) quality_flags = compute_quality_flags(study, comparison, confidence) @@ -52,7 +56,7 @@ def harmonize_metabolomics( quality_flags = sorted(set(quality_flags) | {"identifier_mapping_uncertain"}) rows.append({ - "evidence_id": make_evidence_id(study["study_id"], comparison["comparison_id"], raw_id, analysis_method), + "evidence_id": make_evidence_id(study["study_id"], comparison["comparison_id"], raw_id, method), "study_id": study["study_id"], "comparison_id": comparison["comparison_id"], "simulated": bool(study.get("simulated", False)), @@ -65,7 +69,7 @@ def harmonize_metabolomics( "molecular_direction": direction, "effect_size": effect_size, "effect_size_type": "log2FoldChange", - "standard_error": None, + "standard_error": standard_error, "ci_lower": None, "ci_upper": None, "p_value": float(r["pvalue"]) if pd.notna(r.get("pvalue")) else None, @@ -77,7 +81,7 @@ def harmonize_metabolomics( "phenotype": comparison["phenotype"], "phenotype_direction": comparison["phenotype_direction"], "resilience_classification": comparison.get("resilience_classification"), - "analysis_method": analysis_method, + "analysis_method": method, "mapping_confidence": confidence, "quality_flags": quality_flags, "source_file": source_ref, diff --git a/tests/fixtures/metabolomics_raw/README.md b/tests/fixtures/metabolomics_raw/README.md new file mode 100644 index 0000000..382b88c --- /dev/null +++ b/tests/fixtures/metabolomics_raw/README.md @@ -0,0 +1,16 @@ +# Metabolomics raw-mode fixture (SIMULATED) + +Synthetic, for CI only. Generated with a fixed seed (20261016) to exercise the +targeted-panel options of `workflows/metabolomics` in `raw_reanalysis` mode: + +- two platforms (`LC`, `GC`), each with one internal standard + (`is_internal_standard: true`) and its own per-sample loading, which + `--metabolomics.normalization internal_standard` must remove; +- 0.0001 below-detection placeholders (`--metabolomics.detection_limit 0.001`): + `LC_treatonly` is detected only in treatment and `LC_lod` nowhere, so both + must be reported as not tested, never given a fold change; +- 12 samples in 6 tanks (two families per tank), fitted as 3 vs 3 tank units + with `--metabolomics.replicate_unit tank`. + +Built-in effects (treatment vs control, log2): `LC_up` +2.0, `LC_down` -1.5, +`GC_up` +1.0; the rest are null. diff --git a/tests/fixtures/metabolomics_raw/features.tsv b/tests/fixtures/metabolomics_raw/features.tsv new file mode 100644 index 0000000..e8885c7 --- /dev/null +++ b/tests/fixtures/metabolomics_raw/features.tsv @@ -0,0 +1,14 @@ +feature_id putative_metabolite_name annotation_confidence_level platform is_internal_standard FA-fed-t1 FB-fed-t1 FA-fed-t2 FB-fed-t2 FA-fed-t3 FB-fed-t3 FA-starved-t1 FB-starved-t1 FA-starved-t2 FB-starved-t2 FA-starved-t3 FB-starved-t3 +LC_IS IS-d4 1 LC True 43885.0266 80440.9303 67093.7945 40298.5295 46238.8252 67480.4343 48142.05 52821.5609 53010.4075 41760.5961 44300.9877 109474.9283 +LC_up alanine 1 LC False 11089.3939 24877.379 18716.5127 10386.3421 14083.3126 16837.6314 57496.9815 52670.4284 53489.9972 46541.2837 48594.1556 118891.1204 +LC_down glycine 1 LC False 5891.3451 12484.1608 9374.4616 5341.2766 6397.9503 8811.2911 2582.5243 2099.4937 2223.9817 2076.6679 1845.8413 5051.0709 +LC_null0 serine 1 LC False 2932.353 5049.6953 3968.2056 2916.9544 3308.9209 4285.7926 3497.0162 3608.5314 3443.2045 2956.3322 2872.7604 7013.0896 +LC_null1 proline 1 LC False 5285.1184 11150.9535 7733.3514 5256.8726 6484.436 8179.7812 7059.535 6722.677 6880.2427 5954.1215 6404.7043 15550.8575 +LC_null2 valine 1 LC False 10689.3366 22479.2551 17280.1377 9940.4151 11966.838 16636.0328 13641.7365 13259.4601 13216.4069 11782.3808 13348.0861 27348.325 +LC_treatonly betaine 1 LC False 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 1628.9951 1881.1079 1556.1737 1608.9791 1638.5873 3322.5415 +LC_lod taurine 1 LC False 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 +GC_IS d27-myristic acid 1 GC True 251142.6677 213708.6239 175421.8827 187036.8763 119839.3273 225773.9875 239619.7282 467811.0328 235074.1173 189326.8496 197687.786 418603.8139 +GC_up succinate 1 GC False 35082.228 35363.3624 28371.6898 25387.9837 17182.9244 33787.2403 62202.8408 134951.8409 67504.2825 43932.5721 56022.8724 107531.8974 +GC_null0 citrate 1 GC False 19005.5274 14397.9727 13755.3619 10534.762 6548.6776 15827.1823 19449.3953 29770.3594 18317.0207 11498.9999 11055.3804 24738.3406 +GC_null1 fumarate 1 GC False 8459.2272 6666.8582 6950.028 6258.8621 3701.115 6642.6989 8310.2743 17170.8838 8637.3201 6694.5009 5663.1875 12127.0383 +GC_null2 malate 1 GC False 4247.5338 3570.5364 3430.4191 3043.936 1869.5805 3951.2139 3952.6537 8155.2236 4099.093 3398.4051 2837.181 6295.4736 diff --git a/tests/fixtures/metabolomics_raw/samplesheet.tsv b/tests/fixtures/metabolomics_raw/samplesheet.tsv new file mode 100644 index 0000000..685d682 --- /dev/null +++ b/tests/fixtures/metabolomics_raw/samplesheet.tsv @@ -0,0 +1,13 @@ +sample_id condition family tank +FA-fed-t1 control FA fed_t1 +FB-fed-t1 control FB fed_t1 +FA-fed-t2 control FA fed_t2 +FB-fed-t2 control FB fed_t2 +FA-fed-t3 control FA fed_t3 +FB-fed-t3 control FB fed_t3 +FA-starved-t1 treatment FA starved_t1 +FB-starved-t1 treatment FB starved_t1 +FA-starved-t2 treatment FA starved_t2 +FB-starved-t2 treatment FB starved_t2 +FA-starved-t3 treatment FA starved_t3 +FB-starved-t3 treatment FB starved_t3 diff --git a/tests/test_demo_pipeline.py b/tests/test_demo_pipeline.py index ca10385..dc1216f 100644 --- a/tests/test_demo_pipeline.py +++ b/tests/test_demo_pipeline.py @@ -72,6 +72,39 @@ def test_harmonize_accepts_an_explicit_comparison(tmp_path): assert df["standard_error"].notna().all() +def test_metabolomics_handoff_carries_standard_error_and_method(tmp_path): + """The metabolomics raw mode writes lfcSE and a per-row analysis_method; + both must reach the evidence record. A processed table without them (the + demo's) keeps the old defaults.""" + import pandas as pd + + from harmonize.core import harmonize_processed_table + + out = tmp_path / "GIGAS_GROW06_hypoxia_vs_control_metabolite_features_standardized.tsv" + pd.DataFrame({ + "feature_id": ["GM_alanine", "GM_glycine"], + "putative_metabolite_name": ["alanine", "glycine"], + "annotation_confidence_level": [1, 1], + "metabolite_class": ["amino_acid", "amino_acid"], + "log2FC": [1.98, -1.61], + "pvalue": [3e-10, 5e-9], + "padj": [3e-9, 2e-8], + "lfcSE": [0.032, 0.041], + "analysis_method": ["targeted_MS_limma_moderated_t"] * 2, + }).to_csv(out, sep="\t", index=False) + df = harmonize_processed_table("GIGAS_GROW06", out, date_generated="2026-01-01", + comparison_id="hypoxia_vs_control") + assert df["standard_error"].tolist() == [0.032, 0.041] + assert set(df["analysis_method"]) == {"targeted_MS_limma_moderated_t"} + + legacy = out.with_name("legacy.tsv") + pd.read_csv(out, sep="\t").drop(columns=["lfcSE", "analysis_method"]).to_csv(legacy, sep="\t", index=False) + df = harmonize_processed_table("GIGAS_GROW06", legacy, date_generated="2026-01-01", + comparison_id="hypoxia_vs_control") + assert df["standard_error"].isna().all() + assert set(df["analysis_method"]) == {"untargeted_LCMS"} + + def test_harmonize_rejects_an_unknown_comparison(tmp_path): import pandas as pd import pytest diff --git a/workflows/metabolomics/README.md b/workflows/metabolomics/README.md index a359c31..93c529d 100644 --- a/workflows/metabolomics/README.md +++ b/workflows/metabolomics/README.md @@ -1,11 +1,12 @@ # AREE metabolomics workflow -STATUS: **`processed_results_harmonization` mode is executed, verified, and -covered by current-version CI; the raw feature-table path is not.** It did not -run before 2026-08-28; see `../../docs/first_raw_reanalysis.md` for the defects -that prevented it. No container has been pulled. The raw branch contains real -script logic, but no raw feature-table fixture or real raw study has exercised -it. +STATUS: **both modes execute.** `processed_results_harmonization` is verified +against the bundled demo. `raw_reanalysis` runs natively (`-profile local`) on +the simulated targeted-panel fixture in `tests/fixtures/metabolomics_raw/` +(checked in CI by `scripts/check_metabolomics_raw_smoke.py`) and on the real +targeted panel of `DUPERRET2025_POMS_STARVE`. No container has been pulled. +See `../../docs/first_raw_reanalysis.md` for the defects that kept the +processed path from running before 2026-08-28. This implements Layer 2.D ("Metabolomics") of the AREE reanalysis workflows described in `CLAUDE.md`: metabolite feature table intake, annotation @@ -94,9 +95,41 @@ given file — this mirrors the rnaseq/proteomics workflows and | raw_reanalysis | `--metabolomics.metabolite_annotation_map` (optional) | TSV: `feature_id`, `putative_metabolite_name`, `annotation_confidence_level`. Used to fill in confidence levels not already present in the feature table. | | raw_reanalysis | `--metabolomics.metabolite_class_map` (optional) | TSV: `putative_metabolite_name`, `metabolite_class`[, `pathway`]. See "Metabolite class/pathway crosswalk" below. | | raw_reanalysis | `--metabolomics.qc_cv_threshold` | Coefficient-of-variation threshold for the pooled-QC-sample QC metric (default `0.30`). | +| raw_reanalysis | `--metabolomics.normalization` | `tic` (default), `median` or `internal_standard`; see "Targeted panels" below. | +| raw_reanalysis | `--metabolomics.detection_limit` | Values at or below this become missing (below detection), not intensities. Default `null`. | +| raw_reanalysis | `--metabolomics.replicate_unit` | Sample-sheet column (e.g. `tank`) whose units are the replicates; log2 values are averaged per unit before limma. Default `null`. | +| raw_reanalysis | `--metabolomics.min_detected_per_arm` | Test a feature only if detected in at least this many replicates per arm (default `2`). | +| raw_reanalysis | `--metabolomics.acquisition` | `untargeted` (default) or `targeted`; written to each row's `analysis_method`. | | processed_results_harmonization | `--metabolomics.processed_results` | Already-shaped differential feature table (see target schema below), or a close variant `STANDARDIZE_OUTPUT` can coerce (tolerated header synonyms are listed in `modules/metabolomics/standardize_output.nf`). | | both | `--metabolomics.study_id`, `--metabolomics.comparison_id` | Identifiers threaded into every output filename and the provenance manifest. | +## Targeted panels + +A targeted MRM/SIM panel differs from an untargeted feature table in ways the +default TIC path ignores, so `raw_reanalysis` has options for it: + +- **Platforms.** An optional `platform` column in the feature table (e.g. one + value per LC or GC method) makes every normalization run within a platform, + so panels on different instruments are never scaled against each other. +- **Internal standards.** Rows with `is_internal_standard` true are used by + `--metabolomics.normalization internal_standard` (each sample's mean log2 + standard area, per platform, is subtracted) and are never tested. + `median` centres each sample on its median detected metabolite instead. +- **Detection limit.** Vendor exports often write a placeholder such as + `0.0001` below the limit of detection. `--metabolomics.detection_limit` + turns those into missing values. Features detected in fewer than + `--metabolomics.min_detected_per_arm` replicates of either arm are not + tested, and `__missingness.tsv` (in + `differential_abundance/`) lists every feature with its detection counts + and status, so nothing disappears silently. +- **Replicate units.** `--metabolomics.replicate_unit` averages samples that + share a unit (e.g. two families pooled from one tank) before the fit, the + metabolomics counterpart of the RNA-seq and methylation option. A unit may + not sit in both arms. + +Raw mode also writes `lfcSE` (limma's posterior standard error) and +`analysis_method`, which `aree harmonize` carries into the evidence table. + ## Metabolite class/pathway crosswalk `data/mappings/` currently contains only `gene_id_crosswalk.tsv` (a @@ -124,6 +157,9 @@ tab-separated, with exactly these columns (matching feature_id putative_metabolite_name annotation_confidence_level metabolite_class log2FC pvalue padj ``` +plus, when the input has them (raw mode always does), `lfcSE` and +`analysis_method`. + `annotation_confidence_level` uses a simplified Metabolomics Standards Initiative (MSI)-style 1-4 scale: diff --git a/workflows/metabolomics/assets/report_template.qmd b/workflows/metabolomics/assets/report_template.qmd index a9241c3..8ccad64 100644 --- a/workflows/metabolomics/assets/report_template.qmd +++ b/workflows/metabolomics/assets/report_template.qmd @@ -84,6 +84,8 @@ mode) are embedded in the manifest under `qc_metrics.upstream_step_qc`; in processed-results mode this section is typically empty (see note above). ```{r qc-metrics, results='asis'} +# Base R has %||% only from 4.4. +if (!exists("%||%")) `%||%` <- function(x, y) if (is.null(x)) y else x qc <- manifest$qc_metrics if (!is.null(qc)) { cat(sprintf( @@ -92,6 +94,22 @@ if (!is.null(qc)) { qc$n_annotation_confidence_level_1_or_2 %||% "NA", qc$pct_annotation_confidence_level_1_or_2 %||% "NA" )) + # Raw mode: what normalization and the detection filter did. + for (step in qc$upstream_step_qc) { + if (!is.null(step$normalization_method)) { + cat(sprintf("- Normalization: %s; detection limit %s; %s value(s) below detection; %s internal standard(s) excluded from testing\n", + step$normalization_method, step$detection_limit %||% "none", + step$n_values_below_detection %||% 0, step$n_internal_standards_excluded %||% 0)) + } + if (!is.null(step$n_features_tested)) { + not_tested <- unlist(step$n_features_not_tested) + cat(sprintf("- Tested %s of %s features on %s vs %s replicates%s; not tested: %s (full list in the missingness report)\n", + step$n_features_tested, step$n_features_in, + step$n_replicates_control, step$n_replicates_treatment, + if (is.null(step$replicate_unit)) "" else sprintf(" (units: %s)", step$replicate_unit), + if (length(not_tested)) paste(sprintf("%s %s", not_tested, names(not_tested)), collapse = ", ") else "none")) + } + } } else { cat("No qc_metrics found in manifest.\n") } diff --git a/workflows/metabolomics/main.nf b/workflows/metabolomics/main.nf index 3aac6e2..c834e0a 100644 --- a/workflows/metabolomics/main.nf +++ b/workflows/metabolomics/main.nf @@ -4,7 +4,8 @@ * ------------------------------------ * STATUS: processed_results_harmonization executes end to end against the * bundled demo and is covered by current-version CI. The raw feature-table - * path and all declared containers remain unexecuted; see README.md. + * path runs natively (-profile local) on a CI fixture and on the targeted + * panel of DUPERRET2025_POMS_STARVE; declared containers remain unpulled. * * Implements CLAUDE.md Layer 2.D (Metabolomics): * metabolite feature table intake -> annotation confidence tracking -> @@ -109,14 +110,19 @@ workflow { comparison_id_ch, ANNOTATION_CONFIDENCE.out.annotated_table.map { it[2] }, sample_sheet_ch, - Channel.value(params.metabolomics.qc_cv_threshold) + Channel.value(params.metabolomics.qc_cv_threshold), + Channel.value(params.metabolomics.normalization), + Channel.value(params.metabolomics.detection_limit == null ? 'null' : params.metabolomics.detection_limit) ) DIFFERENTIAL_ABUNDANCE( study_id_ch, comparison_id_ch, NORMALIZE_QC.out.normalized_table.map { it[2] }, - sample_sheet_ch + sample_sheet_ch, + Channel.value(params.metabolomics.replicate_unit ?: ''), + Channel.value(params.metabolomics.min_detected_per_arm), + Channel.value(params.metabolomics.acquisition) ) PATHWAY_MAPPING( @@ -131,6 +137,7 @@ workflow { qc_json_for_manifest = FEATURE_TABLE_INTAKE.out.intake_qc .mix(ANNOTATION_CONFIDENCE.out.annotation_qc) .mix(NORMALIZE_QC.out.normalization_qc) + .mix(DIFFERENTIAL_ABUNDANCE.out.da_qc) .mix(PATHWAY_MAPPING.out.mapping_qc) .collect() diff --git a/workflows/metabolomics/nextflow.config b/workflows/metabolomics/nextflow.config index 59e6d65..828bf80 100644 --- a/workflows/metabolomics/nextflow.config +++ b/workflows/metabolomics/nextflow.config @@ -48,6 +48,31 @@ params { // Coefficient-of-variation threshold for pooled-QC-sample QC metric. qc_cv_threshold = 0.30 + // --- raw_reanalysis: targeted panels and replicate structure --- + // 'tic' (default: scale each sample to the median total, log2(x + 1)), + // 'median' (log2, then centre each sample on its median), or + // 'internal_standard' (log2, then subtract each sample's mean log2 + // internal-standard area). All three run within each value of the + // feature table's optional `platform` column, so panels acquired on + // different instruments are never scaled against each other. Rows with + // `is_internal_standard` true are used only by 'internal_standard' and + // are never tested. + normalization = 'tic' + // Values at or below this are below detection (e.g. a 0.0001 LOD + // placeholder) and become missing, not zero. null keeps the 'tic' + // behaviour of treating them as intensities; 'median' and + // 'internal_standard' always treat non-positive values as missing. + detection_limit = null + // Sample-sheet column (e.g. tank) whose units are the replicates: log2 + // values are averaged per unit before limma. null fits samples. + replicate_unit = null + // A feature is tested only if it is detected in at least this many + // replicates (units, if replicate_unit is set) in each arm. Others are + // listed in the missingness report, not in the results. + min_detected_per_arm = 2 + // 'untargeted' | 'targeted': recorded as the evidence analysis_method. + acquisition = 'untargeted' + // Identifiers threaded through every process/output filename and // into the provenance manifest. Overridden by config/demo.config in // the demo run.