[CW-7172] Improve testing and reporting on differential analysis analysis

README.md

@@ -397,6 +397,12 @@ Output files may be aggregated including information for all samples or provided
 | Differential transcript usage gene summary | de_analysis/{{ contrast }}/results_dtu_gene.tsv | Gene-level DTU summary for one contrast. | aggregated |
 | DEXSeq results | de_analysis/{{ contrast }}/results_dexseq.tsv | Full DEXSeq result table for one contrast. | aggregated |
 | Differential transcript usage plots | de_analysis/{{ contrast }}/results_dtu.pdf | PDF plots generated during DEXSeq analysis for one contrast. | aggregated |
+| Differential analysis QC summary | de_analysis/de_qc_stats.json | Structured DE/DTU QC summary. Use analysis_fallbacks for aggregate counts, and each contrast's deseq2_dispersion_fallback, dexseq_dispersion_method, and dexseq_covariates_dropped fields for interpretation. | aggregated |
+| Differential analysis text summary | de_analysis/de_overall_summary.txt | Human-readable DE/DTU run summary across all contrasts. | aggregated |
+| Per-contrast QC summary | de_analysis/{{ contrast }}/contrast_qc_summary.txt | Human-readable per-contrast DE/DTU QC summary including sample counts and key significance totals. | aggregated |
+| DESeq2 fallback diagnostic | de_analysis/DESeq2_dispersion_fallback_{{ contrast }}.txt | Diagnostic details when DESeq2 falls back to gene-wise dispersion estimation. | aggregated |
+| DTU failure diagnostic | de_analysis/{{ contrast }}/DTU_ANALYSIS_FAILED.txt | Diagnostic details when DEXSeq fails for a contrast. | aggregated |
+| Multiple-testing warning | de_analysis/MULTIPLE_TESTING_WARNING.txt | Family-wise error-rate note generated when multiple contrasts are tested. | aggregated |
 | IGV configuration | igv.json | JSON configuration for viewing the aligned BAMs in IGV. | aggregated |
 | Reference FASTA index | igv_reference/{{ ref_genome_file }}.fai | FAI index for the reference genome published for IGV. | aggregated |
 | Reference GZI index | igv_reference/{{ ref_genome_file }}.gzi | GZI index for a compressed reference genome published for IGV. | aggregated |

bin/workflow_glue/report.py

@@ -171,6 +171,64 @@ def _create_warning_banner(message, level="warning"):
         raw(message)
 
 
+def _as_string_list(value):
+    """Normalize optional values to a compact list of strings."""
+    if value is None or value == "none":
+        return []
+    if isinstance(value, (list, tuple, set)):
+        return [str(item) for item in value if item not in (None, "")]
+    if isinstance(value, str):
+        return [value] if value else []
+    return [str(value)]
+
+
+def _collect_de_method_rows(de_qc):
+    """Build per-contrast method rows and warning metadata."""
+    rows = []
+    deseq2_gene_wise = []
+    dexseq_gene_wise = []
+    dexseq_covariate_drops = []
+
+    for contrast_name, contrast_data in de_qc.get("contrasts", {}).items():
+        fallback = contrast_data.get("deseq2_dispersion_fallback") or {}
+        fallback_applied = bool(fallback.get("applied", False))
+        deseq2_method = fallback.get("method_used")
+
+        if not deseq2_method:
+            deseq2_method = "gene-wise" if fallback_applied else "parametric"
+
+        if deseq2_method == "gene-wise":
+            deseq2_gene_wise.append(contrast_name)
+
+        dexseq_method = contrast_data.get("dexseq_dispersion_method") or "parametric"
+        if dexseq_method == "gene-wise":
+            dexseq_gene_wise.append(contrast_name)
+
+        dropped_covariates = _as_string_list(
+            contrast_data.get("dexseq_covariates_dropped")
+        )
+        if dropped_covariates:
+            dexseq_covariate_drops.append((contrast_name, dropped_covariates))
+
+        rows.append(
+            {
+                "Contrast": contrast_name,
+                "DESeq2 dispersion": (
+                    f"{deseq2_method} (fallback)"
+                    if fallback_applied
+                    else deseq2_method
+                ),
+                "DEXSeq dispersion": dexseq_method,
+                "DEXSeq covariates dropped": (
+                    ", ".join(dropped_covariates) if dropped_covariates else "none"
+                ),
+                "DTU status": contrast_data.get("dtu_status", "N/A"),
+            }
+        )
+
+    return rows, deseq2_gene_wise, dexseq_gene_wise, dexseq_covariate_drops
+
+
 def main(args):
     """Run the report entry point."""
     logger = get_named_logger("Report")
@@ -471,13 +529,21 @@ def main(args):
             ):
                 # Check for critical warnings
                 has_warnings = False
-
+                sample_size_warnings = _as_string_list(
-                if (
                     de_qc.get("sample_size_warnings")
-                    and de_qc["sample_size_warnings"] != "none"
+                )
-                ):
+                (
+                    method_rows,
+                    deseq2_gene_wise,
+                    dexseq_gene_wise,
+                    dexseq_covariate_drops,
+                ) = _collect_de_method_rows(de_qc)
+
+                if sample_size_warnings:
                     _create_warning_banner(
-                        f"Sample Size Warning: {de_qc['sample_size_warnings']}. "
+                        "Sample Size Warning: "
+                        + "; ".join(sample_size_warnings)
+                        + ". "
                         "Underpowered designs may have reduced statistical "
                         "power and increased false negative rate.",
                         level="warning",
@@ -491,24 +557,31 @@ def main(args):
                         level="info",
                     )
 
-                # Check for dispersion fallbacks
+                if deseq2_gene_wise or dexseq_gene_wise:
-                dispersion_fallbacks = []
+                    gene_wise_details = []
-                for contrast_name, contrast_data in de_qc.get(
+                    if deseq2_gene_wise:
-                    "contrasts", {}
+                        gene_wise_details.append(
-                ).items():
+                            "DESeq2: " + ", ".join(sorted(deseq2_gene_wise))
-                    dispersion_file = (
+                        )
-                        Path(args.de_dir)
+                    if dexseq_gene_wise:
-                        / f"DESeq2_dispersion_fallback_{contrast_name}.txt"
+                        gene_wise_details.append(
+                            "DEXSeq: " + ", ".join(sorted(dexseq_gene_wise))
                         )
-                    if dispersion_file.exists():
-                        dispersion_fallbacks.append(contrast_name)
-
-                if dispersion_fallbacks:
                     _create_warning_banner(
-                        "Dispersion Estimation Fallback: "
+                        "Gene-wise dispersion fallback used (reduced power). "
-                        f"{len(dispersion_fallbacks)} contrast(s) used "
+                        + " ".join(gene_wise_details),
-                        "gene-wise dispersion (reduced power). "
+                        level="warning",
-                        f"Affected: {', '.join(dispersion_fallbacks)}",
+                    )
+                    has_warnings = True
+
+                if dexseq_covariate_drops:
+                    drop_details = [
+                        f"{contrast} ({', '.join(columns)})"
+                        for contrast, columns in sorted(dexseq_covariate_drops)
+                    ]
+                    _create_warning_banner(
+                        "DEXSeq covariates dropped due to rank-deficient design. "
+                        f"Affected: {'; '.join(drop_details)}",
                         level="warning",
                     )
                     has_warnings = True
@@ -534,12 +607,8 @@ def main(args):
 
                 # Experimental design summary
                 with h3("Experimental Design"):
-                    covariates = de_qc.get("covariates", [])
+                    covariates = _as_string_list(de_qc.get("covariates"))
-                    covariates_value = (
+                    covariates_value = ", ".join(covariates) if covariates else "none"
-                        ", ".join(covariates)
-                        if de_qc.get("covariates") != "none"
-                        else "none"
-                    )
                     design_stats = pd.DataFrame(
                         [
                             ("Total samples", de_qc.get("total_samples", 0)),
@@ -591,6 +660,17 @@ def main(args):
                             use_index=False,
                         )
 
+                with h3("Statistical Methods & Warnings"):
+                    if method_rows:
+                        method_df = pd.DataFrame(method_rows)
+                        DataTable.from_pandas(
+                            method_df,
+                            paging=False,
+                            use_index=False,
+                        )
+                    else:
+                        p("No contrast-level QC metadata was found.")
+
                 # Per-contrast summary
                 if "contrasts" in de_qc:
                     with h3("Results Summary by Contrast"):
@@ -637,22 +717,35 @@ def main(args):
                 if has_warnings:
                     with h3("Quality Warnings Summary"):
                         warnings_data = []
-                        if (
+                        if sample_size_warnings:
-                            de_qc.get("sample_size_warnings")
-                            and de_qc["sample_size_warnings"] != "none"
-                        ):
                             warnings_data.append(
                                 {
                                     "Warning Type": "Sample Size",
-                                    "Details": de_qc["sample_size_warnings"],
+                                    "Details": "; ".join(sample_size_warnings),
                                 }
                             )
-                        if dispersion_fallbacks:
+                        if deseq2_gene_wise or dexseq_gene_wise:
+                            engines = []
+                            if deseq2_gene_wise:
+                                engines.append(
+                                    f"DESeq2 ({len(deseq2_gene_wise)} contrasts)"
+                                )
+                            if dexseq_gene_wise:
+                                engines.append(
+                                    f"DEXSeq ({len(dexseq_gene_wise)} contrasts)"
+                                )
                             warnings_data.append(
                                 {
-                                    "Warning Type": "Dispersion Estimation",
+                                    "Warning Type": "Gene-wise Dispersion Fallback",
+                                    "Details": "; ".join(engines),
+                                }
+                            )
+                        if dexseq_covariate_drops:
+                            warnings_data.append(
+                                {
+                                    "Warning Type": "DEXSeq Covariates Dropped",
                                     "Details": (
-                                        f"{len(dispersion_fallbacks)} "
+                                        f"{len(dexseq_covariate_drops)} "
                                         "contrasts affected"
                                     ),
                                 }

bin/workflow_glue/tests/test_report.py

@@ -42,6 +42,90 @@ def _write(path, text):
     return path
 
 
+def _build_report_args(tmp_path, de_qc=None):
+    """Create minimal report inputs, optionally including DE QC JSON."""
+    metadata = _write(
+        tmp_path / "metadata.json",
+        json.dumps([{"alias": "sampleA", "has_stats": False}]),
+    )
+    params = _write(tmp_path / "params.json", "{}")
+    versions = tmp_path / "versions"
+    versions.mkdir()
+    _write(versions / "versions.txt", "tool,1.0\n")
+
+    cohort = tmp_path / "cohort"
+    cohort.mkdir()
+    reference = cohort / "reference"
+    reference.mkdir()
+    _write(
+        reference / "annotation_reference_summary.json",
+        json.dumps(
+            {
+                "seqname_overlap": ["chr1"],
+                "only_in_annotation": [],
+                "only_in_reference": [],
+                "annotation": {
+                    "kept_records": 10,
+                    "excluded_unstranded_records": 0,
+                    "sanitised_attribute_records": 0,
+                },
+                "warnings": [],
+            }
+        ),
+    )
+
+    samples = tmp_path / "samples"
+    samples.mkdir()
+    sqanti = tmp_path / "sqanti"
+    sqanti.mkdir()
+    alignment_stats = tmp_path / "alignment_stats"
+    alignment_stats.mkdir()
+    (samples / "OPTIONAL_FILE").touch()
+    (sqanti / "OPTIONAL_FILE").touch()
+    (alignment_stats / "OPTIONAL_FILE").touch()
+
+    de_dir = None
+    if de_qc is not None:
+        de_dir = tmp_path / "de_analysis"
+        de_dir.mkdir()
+        _write(de_dir / "de_qc_stats.json", json.dumps(de_qc))
+        for contrast_name in de_qc.get("contrasts", {}):
+            contrast_dir = de_dir / contrast_name
+            contrast_dir.mkdir()
+            _write(
+                contrast_dir / "results_dge.tsv",
+                "GENEID\tlog2FoldChange\tpadj\n"
+                "gene1\t1.0\t0.01\n",
+            )
+            _write(
+                contrast_dir / "results_dtu_transcript.tsv",
+                "featureID\tgroupID\tpadj\n"
+                "tx1\tgene1\t0.05\n",
+            )
+
+    out_report = tmp_path / "wf-transcriptomes-report.html"
+    argv = [
+        str(out_report),
+        "--metadata",
+        str(metadata),
+        "--alignment_stats_dir",
+        str(alignment_stats),
+        "--cohort_dir",
+        str(cohort),
+        "--samples_dir",
+        str(samples),
+        "--sqanti_dir",
+        str(sqanti),
+        "--versions",
+        str(versions),
+        "--params",
+        str(params),
+    ]
+    if de_dir is not None:
+        argv.extend(["--de_dir", str(de_dir)])
+    return report.argparser().parse_args(argv), out_report
+
+
 def test_report_main_accepts_optional_file_sentinels(monkeypatch, tmp_path):
     """The report entry point should tolerate null-object sentinel files."""
     tables = []
@@ -150,3 +234,153 @@ def test_pychopper_tables_uses_sample_directory_names(tmp_path):
         "Full length",
         "Unclassified",
     ]
+
+
+def test_report_main_renders_statistical_methods_and_warnings(
+    monkeypatch,
+    tmp_path,
+):
+    """DE/DTU QC report renders fallback methods and warning banners."""
+    tables = []
+    headings = []
+    banners = []
+
+    monkeypatch.setattr(report.labs, "LabsReport", _FakeReport)
+    monkeypatch.setattr(report, "Tabs", _FakeTabs)
+    monkeypatch.setattr(report, "p", lambda *args, **kwargs: None)
+    monkeypatch.setattr(report, "pre", lambda *args, **kwargs: None)
+    monkeypatch.setattr(report.fastcat, "SeqSummary", lambda *args, **kwargs: None)
+    monkeypatch.setattr(
+        report,
+        "h3",
+        lambda label: (headings.append(label), _NullContext())[1],
+    )
+    monkeypatch.setattr(
+        report,
+        "_create_warning_banner",
+        lambda message, level="warning": banners.append((level, message)),
+    )
+    monkeypatch.setattr(
+        report.DataTable,
+        "from_pandas",
+        staticmethod(lambda table, *args, **kwargs: tables.append(table.copy())),
+    )
+
+    de_qc = {
+        "total_samples": 6,
+        "condition_column": "condition",
+        "reference_level": "control",
+        "covariates": ["batch"],
+        "num_contrasts": 2,
+        "sample_size_warnings": "none",
+        "samples_per_group": {"control": 3, "treated": 3},
+        "contrasts": {
+            "condition_treated_vs_control": {
+                "n_target": 3,
+                "n_reference": 3,
+                "dge_significant_fdr05": 10,
+                "dge_upregulated": 6,
+                "dge_downregulated": 4,
+                "dtu_status": "SUCCESS",
+                "dtu_significant_genes": 2,
+                "deseq2_dispersion_fallback": {
+                    "applied": True,
+                    "method_used": "gene-wise",
+                    "reason": "recoverable",
+                    "diagnostic_file": (
+                        "DESeq2_dispersion_fallback_"
+                        "condition_treated_vs_control.txt"
+                    ),
+                },
+                "dexseq_dispersion_method": "local",
+                "dexseq_covariates_dropped": ["batch"],
+            },
+            "condition_treated2_vs_control": {
+                "n_target": 3,
+                "n_reference": 3,
+                "dge_significant_fdr05": 4,
+                "dge_upregulated": 3,
+                "dge_downregulated": 1,
+                "dtu_status": "FAILED",
+            },
+        },
+    }
+
+    args, out_report = _build_report_args(tmp_path, de_qc=de_qc)
+    report.main(args)
+
+    assert out_report.exists()
+    assert "Statistical Methods & Warnings" in headings
+    assert any("DESeq2 dispersion" in table.columns for table in tables)
+    assert any(
+        "gene-wise (fallback)" in table.to_string()
+        for table in tables
+        if "DESeq2 dispersion" in table.columns
+    )
+    assert any(
+        "batch" in table.to_string()
+        for table in tables
+        if "DEXSeq covariates dropped" in table.columns
+    )
+    assert any("gene-wise dispersion fallback" in msg.lower() for _, msg in banners)
+    assert any("covariates dropped" in msg.lower() for _, msg in banners)
+    assert any(
+        level == "danger" and "DTU Analysis Failed" in msg
+        for level, msg in banners
+    )
+
+
+def test_report_main_tolerates_missing_statistical_fields(monkeypatch, tmp_path):
+    """Older DE QC JSON without new fallback fields should still render."""
+    tables = []
+    headings = []
+
+    monkeypatch.setattr(report.labs, "LabsReport", _FakeReport)
+    monkeypatch.setattr(report, "Tabs", _FakeTabs)
+    monkeypatch.setattr(report, "p", lambda *args, **kwargs: None)
+    monkeypatch.setattr(report, "pre", lambda *args, **kwargs: None)
+    monkeypatch.setattr(report.fastcat, "SeqSummary", lambda *args, **kwargs: None)
+    monkeypatch.setattr(
+        report,
+        "h3",
+        lambda label: (headings.append(label), _NullContext())[1],
+    )
+    monkeypatch.setattr(report, "_create_warning_banner", lambda *args, **kwargs: None)
+    monkeypatch.setattr(
+        report.DataTable,
+        "from_pandas",
+        staticmethod(lambda table, *args, **kwargs: tables.append(table.copy())),
+    )
+
+    legacy_de_qc = {
+        "total_samples": 4,
+        "condition_column": "condition",
+        "reference_level": "control",
+        "covariates": "none",
+        "num_contrasts": 1,
+        "sample_size_warnings": "none",
+        "samples_per_group": {"control": 2, "treated": 2},
+        "contrasts": {
+            "condition_treated_vs_control": {
+                "n_target": 2,
+                "n_reference": 2,
+                "dge_significant_fdr05": 1,
+                "dge_upregulated": 1,
+                "dge_downregulated": 0,
+                "dtu_status": "SUCCESS",
+            }
+        },
+    }
+
+    args, out_report = _build_report_args(tmp_path, de_qc=legacy_de_qc)
+    report.main(args)
+
+    assert out_report.exists()
+    assert "Statistical Methods & Warnings" in headings
+    method_tables = [
+        table
+        for table in tables
+        if "DESeq2 dispersion" in table.columns
+    ]
+    assert method_tables
+    assert "parametric" in method_tables[0].to_string()

bin/workflow_glue_r/R/de_analysis.R

@@ -145,8 +145,23 @@ de_set_dispersions <- function(object, value) {
     setter(object, value = value)
 }
 
-de_run_deseq_with_fallback <- function(dds, contrast_name, out_dir) {
+de_extract_disp_gene_est <- function(object) {
-    tryCatch(
+    S4Vectors::mcols(object)$dispGeneEst
+}
+
+de_run_deseq_with_fallback <- function(
+    dds,
+    contrast_name,
+    out_dir
+) {
+    fallback_info <- list(
+        applied = FALSE,
+        method_used = "parametric",
+        reason = NULL,
+        diagnostic_file = NULL
+    )
+
+    de_out <- tryCatch(
         DESeq2::DESeq(dds, quiet = TRUE),
         error = function(err) {
             if (!grepl(
@@ -171,7 +186,19 @@ de_run_deseq_with_fallback <- function(dds, contrast_name, out_dir) {
 
             dds <- DESeq2::estimateSizeFactors(dds)
             dds <- DESeq2::estimateDispersionsGeneEst(dds)
-            dds <- de_set_dispersions(dds, S4Vectors::mcols(dds)$dispGeneEst)
+            dds <- de_set_dispersions(dds, de_extract_disp_gene_est(dds))
+
+            dispersion_values <- suppressWarnings(as.numeric(DESeq2::dispersions(dds)))
+            dispersion_values <- dispersion_values[is.finite(dispersion_values)]
+            dispersion_range <- if (length(dispersion_values) > 0) {
+                sprintf(
+                    "Dispersion range: %.3f to %.3f",
+                    min(dispersion_values),
+                    max(dispersion_values)
+                )
+            } else {
+                "Dispersion range: unavailable"
+            }
 
             diag_content <- c(
                 "DESeq2 Dispersion Estimation Fallback Applied",
@@ -181,11 +208,7 @@ de_run_deseq_with_fallback <- function(dds, contrast_name, out_dir) {
                 sprintf("Contrast: %s", contrast_name),
                 sprintf("Samples: %d", ncol(dds)),
                 sprintf("Genes tested: %d", nrow(dds)),
-                sprintf(
+                dispersion_range,
-                    "Dispersion range: %.3f to %.3f",
-                    min(DESeq2::dispersions(dds)),
-                    max(DESeq2::dispersions(dds))
-                ),
                 "",
                 "WHAT HAPPENED:",
                 "  Curve fitting failed. Using gene-wise dispersion estimates.",
@@ -215,15 +238,31 @@ de_run_deseq_with_fallback <- function(dds, contrast_name, out_dir) {
                 )
             )
             writeLines(diag_content, diag_file)
+            fallback_info <<- list(
+                applied = TRUE,
+                method_used = "gene-wise",
+                reason = conditionMessage(err),
+                diagnostic_file = basename(diag_file)
+            )
 
             DESeq2::nbinomWaldTest(dds)
         }
     )
+    list(dds = de_out, deseq2_dispersion_fallback = fallback_info)
 }
 
-de_estimate_dispersions_with_fallback <- function(object, context_label, allow_gene_est = TRUE) {
+de_estimate_dispersions_with_fallback <- function(
+    object,
+    context_label,
+    allow_gene_est = TRUE
+) {
     tryCatch(
-        DESeq2::estimateDispersions(object),
+        list(
+            object = DESeq2::estimateDispersions(object),
+            method_used = "parametric",
+            fallback_applied = FALSE,
+            reason = NULL
+        ),
         error = function(err) {
             if (!grepl(
                 "all gene-wise dispersion estimates are within 2 orders of magnitude",
@@ -233,9 +272,18 @@ de_estimate_dispersions_with_fallback <- function(object, context_label, allow_g
                 stop(err)
             }
 
-            message(context_label, " dispersion fitting failed; retrying with fitType='local'.")
+            primary_reason <- conditionMessage(err)
+            message(
+                context_label,
+                " dispersion fitting failed; retrying with fitType='local'."
+            )
             tryCatch(
-                DESeq2::estimateDispersions(object, fitType = "local"),
+                list(
+                    object = DESeq2::estimateDispersions(object, fitType = "local"),
+                    method_used = "local",
+                    fallback_applied = TRUE,
+                    reason = primary_reason
+                ),
                 error = function(local_err) {
                     if (!grepl(
                         "all gene-wise dispersion estimates are within 2 orders of magnitude",
@@ -245,9 +293,17 @@ de_estimate_dispersions_with_fallback <- function(object, context_label, allow_g
                         stop(local_err)
                     }
 
-                    message(context_label, " local-fit dispersion retry failed; retrying with fitType='mean'.")
+                    message(
+                        context_label,
+                        " local-fit dispersion retry failed; retrying with fitType='mean'."
+                    )
                     tryCatch(
-                        DESeq2::estimateDispersions(object, fitType = "mean"),
+                        list(
+                            object = DESeq2::estimateDispersions(object, fitType = "mean"),
+                            method_used = "mean",
+                            fallback_applied = TRUE,
+                            reason = primary_reason
+                        ),
                         error = function(mean_err) {
                             if (!grepl(
                                 "all gene-wise dispersion estimates are within 2 orders of magnitude",
@@ -265,8 +321,13 @@ de_estimate_dispersions_with_fallback <- function(object, context_label, allow_g
                                 " mean-fit dispersion retry failed; falling back to gene-wise dispersion estimates."
                             )
                             object <- DESeq2::estimateDispersionsGeneEst(object)
-                            object <- de_set_dispersions(object, S4Vectors::mcols(object)$dispGeneEst)
+                            object <- de_set_dispersions(object, de_extract_disp_gene_est(object))
-                            object
+                            list(
+                                object = object,
+                                method_used = "gene-wise",
+                                fallback_applied = TRUE,
+                                reason = primary_reason
+                            )
                         }
                     )
                 }
@@ -315,13 +376,19 @@ de_run_deseq2_result <- function(
         colData = coldata,
         design = design_formula
     )
-    dds <- de_run_deseq_with_fallback(dds, contrast_name, out_dir)
+    deseq_run <- de_run_deseq_with_fallback(dds, contrast_name, out_dir)
+    dds <- deseq_run$dds
+    deseq2_dispersion_fallback <- deseq_run$deseq2_dispersion_fallback
     result <- DESeq2::results(
         dds,
         contrast = c(condition_column, target_level, reference_level),
         independentFiltering = TRUE
     )
-    list(dds = dds, result = result)
+    list(
+        dds = dds,
+        result = result,
+        deseq2_dispersion_fallback = deseq2_dispersion_fallback
+    )
 }
 
 de_run_dexseq_result <- function(
@@ -336,6 +403,7 @@ de_run_dexseq_result <- function(
     for (covariate in covariates) {
         coldata[[covariate]] <- factor(coldata[[covariate]])
     }
+    dropped_covariates <- character(0)
 
     run_inner <- function(active_covariates) {
         covariate_exon_terms <- if (length(active_covariates) > 0) {
@@ -357,11 +425,29 @@ de_run_dexseq_result <- function(
                 groupID = tx_meta$GENEID
             )
             dxd <- DESeq2::estimateSizeFactors(dxd)
-            dxd <- de_estimate_dispersions_with_fallback(dxd, "DEXSeq", allow_gene_est = TRUE)
+            dispersion_result <- de_estimate_dispersions_with_fallback(
+                dxd,
+                "DEXSeq",
+                allow_gene_est = TRUE
+            )
+            if (is.list(dispersion_result) && !is.null(dispersion_result$object)) {
+                dxd <- dispersion_result$object
+                dispersion_method <- dispersion_result$method_used
+                dispersion_reason <- dispersion_result$reason
+            } else {
+                dxd <- dispersion_result
+                dispersion_method <- "parametric"
+                dispersion_reason <- NULL
+            }
             dxd <- DEXSeq::testForDEU(dxd, reducedModel = reduced_formula)
             dxd <- DEXSeq::estimateExonFoldChanges(dxd, fitExpToVar = condition_column)
             dxr <- DEXSeq::DEXSeqResults(dxd, independentFiltering = FALSE)
-            list(dxd = dxd, dxr = dxr)
+            list(
+                dxd = dxd,
+                dxr = dxr,
+                dexseq_dispersion_method = dispersion_method,
+                dexseq_dispersion_reason = dispersion_reason
+            )
         }, error = function(err) {
             if (length(active_covariates) == 0 || !grepl(
                 "model matrix is not full rank",
@@ -373,6 +459,7 @@ de_run_dexseq_result <- function(
 
             dropped_covariate <- tail(active_covariates, 1)
             kept_covariates <- head(active_covariates, -1)
+            dropped_covariates <<- c(dropped_covariates, dropped_covariate)
             message(
                 "DEXSeq design was not full rank with covariate '",
                 dropped_covariate,
@@ -382,20 +469,12 @@ de_run_dexseq_result <- function(
         })
     }
 
-    run_inner(covariates)
+    result <- run_inner(covariates)
+    result$dexseq_covariates_dropped <- dropped_covariates
+    result
 }
 
-main_run_de_analysis <- function(
+main_run_de_analysis <- function(argv) {
-    argv,
-    deseq_runner = de_run_deseq2_result,
-    dexseq_runner = de_run_dexseq_result,
-    pdf_fn = grDevices::pdf,
-    dev_off_fn = grDevices::dev.off,
-    plot_ma_fn = DESeq2::plotMA,
-    plot_disp_fn = DESeq2::plotDispEsts,
-    per_gene_q_fn = DEXSeq::perGeneQValue,
-    placeholder_pdf_fn = de_write_placeholder_pdf
-) {
     set.seed(42)
     dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE)
 
@@ -512,7 +591,15 @@ main_run_de_analysis <- function(
             reference_level = reference_level,
             n_samples = nrow(contrast_samples),
             n_target = sum(contrast_samples[[argv$condition_column]] == target_level),
-            n_reference = sum(contrast_samples[[argv$condition_column]] == reference_level)
+            n_reference = sum(contrast_samples[[argv$condition_column]] == reference_level),
+            deseq2_dispersion_fallback = list(
+                applied = FALSE,
+                method_used = "parametric",
+                reason = NULL,
+                diagnostic_file = NULL
+            ),
+            dexseq_dispersion_method = "parametric",
+            dexseq_covariates_dropped = list()
         )
 
         if (nrow(contrast_samples) < 6) {
@@ -528,7 +615,7 @@ main_run_de_analysis <- function(
         contrast_qc$genes_tested <- nrow(gene_counts)
         contrast_qc$transcripts_tested <- nrow(tx_counts)
 
-        dge_run <- deseq_runner(
+        dge_run <- de_run_deseq2_result(
             gene_counts,
             contrast_samples,
             target_level,
@@ -538,6 +625,20 @@ main_run_de_analysis <- function(
             argv$out_dir,
             contrast_name
         )
+        if (!is.null(dge_run$deseq2_dispersion_fallback)) {
+            fallback <- dge_run$deseq2_dispersion_fallback
+            fallback_applied <- isTRUE(fallback$applied)
+            fallback_method <- fallback$method_used
+            if (is.null(fallback_method) || identical(fallback_method, "")) {
+                fallback_method <- if (fallback_applied) "gene-wise" else "parametric"
+            }
+            contrast_qc$deseq2_dispersion_fallback <- list(
+                applied = fallback_applied,
+                method_used = fallback_method,
+                reason = fallback$reason,
+                diagnostic_file = fallback$diagnostic_file
+            )
+        }
         dge_res <- as.data.frame(dge_run$result)
         dge_res$GENEID <- rownames(dge_res)
         dge_res <- merge(gene_meta, dge_res, by = "GENEID", all.y = TRUE, sort = FALSE)
@@ -563,12 +664,12 @@ main_run_de_analysis <- function(
             row.names = FALSE
         )
 
-        pdf_fn(file.path(contrast_dir, "results_dge.pdf"))
+        grDevices::pdf(file.path(contrast_dir, "results_dge.pdf"))
-        plot_ma_fn(dge_run$result)
+        DESeq2::plotMA(dge_run$result)
-        dev_off_fn()
+        grDevices::dev.off()
 
         dex_res <- tryCatch(
-            dexseq_runner(
+            de_run_dexseq_result(
                 tx_counts,
                 tx_meta,
                 contrast_samples,
@@ -628,7 +729,7 @@ main_run_de_analysis <- function(
             ))
             tx_dtu <- dex_df
             gene_dtu <- workflow_glue_r_empty_tsv(c("GENEID", "qval"))
-            placeholder_pdf_fn(
+            de_write_placeholder_pdf(
                 file.path(contrast_dir, "results_dtu.pdf"),
                 "DEXSeq did not converge for this contrast.\nSee DTU_ANALYSIS_FAILED.txt for details."
             )
@@ -636,6 +737,12 @@ main_run_de_analysis <- function(
             contrast_qc$dtu_significant_transcripts <- 0
             contrast_qc$dtu_significant_genes <- 0
         } else {
+            if (!is.null(dex_res$dexseq_dispersion_method)) {
+                contrast_qc$dexseq_dispersion_method <- dex_res$dexseq_dispersion_method
+            }
+            if (!is.null(dex_res$dexseq_covariates_dropped)) {
+                contrast_qc$dexseq_covariates_dropped <- as.list(dex_res$dexseq_covariates_dropped)
+            }
             dex_df <- as.data.frame(dex_res$dxr)
             dex_df <- workflow_glue_r_normalise_tsv_df(dex_df)
             tx_dtu <- dex_df[, intersect(
@@ -644,7 +751,7 @@ main_run_de_analysis <- function(
             ), drop = FALSE]
             tx_dtu <- workflow_glue_r_normalise_tsv_df(tx_dtu)
 
-            gene_q <- per_gene_q_fn(dex_res$dxr)
+            gene_q <- DEXSeq::perGeneQValue(dex_res$dxr)
             gene_dtu <- data.frame(
                 GENEID = names(gene_q),
                 qval = unname(gene_q),
@@ -655,10 +762,10 @@ main_run_de_analysis <- function(
             contrast_qc$dtu_significant_transcripts <- sum(tx_dtu$padj < 0.05, na.rm = TRUE)
             contrast_qc$dtu_significant_genes <- sum(gene_dtu$qval < 0.05, na.rm = TRUE)
 
-            pdf_fn(file.path(contrast_dir, "results_dtu.pdf"))
+            grDevices::pdf(file.path(contrast_dir, "results_dtu.pdf"))
-            plot_ma_fn(dex_res$dxr, cex = 0.8, alpha = 0.05)
+            DESeq2::plotMA(dex_res$dxr, cex = 0.8, alpha = 0.05)
-            plot_disp_fn(dex_res$dxd)
+            DESeq2::plotDispEsts(dex_res$dxd)
-            dev_off_fn()
+            grDevices::dev.off()
         }
 
         utils::write.table(
@@ -727,6 +834,48 @@ main_run_de_analysis <- function(
         de_qc_stats$contrasts[[contrast_name]] <- contrast_qc
     }
 
+    deseq2_dispersion_fallbacks <- names(Filter(
+        function(cqc) isTRUE(cqc$deseq2_dispersion_fallback$applied),
+        de_qc_stats$contrasts
+    ))
+    deseq2_gene_wise <- names(Filter(
+        function(cqc) identical(cqc$deseq2_dispersion_fallback$method_used, "gene-wise"),
+        de_qc_stats$contrasts
+    ))
+    dexseq_non_parametric <- names(Filter(
+        function(cqc) {
+            method <- cqc$dexseq_dispersion_method
+            !is.null(method) && !identical(method, "parametric")
+        },
+        de_qc_stats$contrasts
+    ))
+    dexseq_gene_wise <- names(Filter(
+        function(cqc) identical(cqc$dexseq_dispersion_method, "gene-wise"),
+        de_qc_stats$contrasts
+    ))
+    dexseq_covariate_drop <- names(Filter(
+        function(cqc) length(cqc$dexseq_covariates_dropped) > 0,
+        de_qc_stats$contrasts
+    ))
+    total_covariates_dropped <- sum(vapply(
+        de_qc_stats$contrasts,
+        function(cqc) length(cqc$dexseq_covariates_dropped),
+        integer(1)
+    ))
+    de_qc_stats$analysis_fallbacks <- list(
+        deseq2_dispersion_fallback_contrasts = length(deseq2_dispersion_fallbacks),
+        deseq2_dispersion_fallback_contrast_names = as.list(deseq2_dispersion_fallbacks),
+        deseq2_gene_wise_contrasts = length(deseq2_gene_wise),
+        deseq2_gene_wise_contrast_names = as.list(deseq2_gene_wise),
+        dexseq_non_parametric_dispersion_contrasts = length(dexseq_non_parametric),
+        dexseq_non_parametric_dispersion_contrast_names = as.list(dexseq_non_parametric),
+        dexseq_gene_wise_dispersion_contrasts = length(dexseq_gene_wise),
+        dexseq_gene_wise_dispersion_contrast_names = as.list(dexseq_gene_wise),
+        dexseq_covariate_drop_contrasts = length(dexseq_covariate_drop),
+        dexseq_covariate_drop_contrast_names = as.list(dexseq_covariate_drop),
+        total_covariates_dropped = total_covariates_dropped
+    )
+
     jsonlite::write_json(
         de_qc_stats,
         file.path(argv$out_dir, "de_qc_stats.json"),

bin/workflow_glue_r/tests/testthat/test_de_analysis.R

@@ -300,11 +300,7 @@ testthat::test_that("transcript SE without GENEID rejected", {
     )
 
     testthat::expect_error(
-        main_run_de_analysis(
+        main_run_de_analysis(argv),
-            argv,
-            deseq_runner = function(...) stop("runner should not be called"),
-            dexseq_runner = function(...) stop("runner should not be called")
-        ),
         "Transcript rowData must contain GENEID"
     )
 })
@@ -347,11 +343,7 @@ testthat::test_that("underspecified designs rejected", {
     )
 
     testthat::expect_error(
-        suppressWarnings(main_run_de_analysis(
+        suppressWarnings(main_run_de_analysis(argv)),
-            argv,
-            deseq_runner = function(...) stop("runner should not be called"),
-            dexseq_runner = function(...) stop("runner should not be called")
-        )),
         "fewer than 2 replicates"
     )
 
@@ -376,11 +368,7 @@ testthat::test_that("underspecified designs rejected", {
     argv$sample_sheet <- single_condition_sheet
     argv$out_dir <- file.path(fixture_dir, "single-out")
     testthat::expect_error(
-        main_run_de_analysis(
+        main_run_de_analysis(argv),
-            argv,
-            deseq_runner = function(...) stop("runner should not be called"),
-            dexseq_runner = function(...) stop("runner should not be called")
-        ),
         "requires at least two condition levels"
     )
 
@@ -390,18 +378,178 @@ testthat::test_that("underspecified designs rejected", {
     argv$sample_sheet <- sample_sheet
     argv$out_dir <- file.path(fixture_dir, "malformed-out")
     testthat::expect_error(
-        main_run_de_analysis(
+        main_run_de_analysis(argv)
-            argv,
+    )
-            deseq_runner = function(...) stop("runner should not be called"),
+})
-            dexseq_runner = function(...) stop("runner should not be called")
+
+###
+# Fallback helpers and metadata wiring
+#
+# Fixture-driven tests for DESeq2/DEXSeq helper behaviour and metadata.
+# These avoid dependency injection and exercise the real package code paths.
+
+testthat::test_that("de_run_deseq_with_fallback returns structured metadata", {
+    testthat::skip_if_not_installed("DESeq2")
+
+    gene_se <- make_test_gene_se()
+    sample_df <- data.frame(
+        alias = colnames(gene_se),
+        condition = rep(c("control", "treated"), each = 3),
+        batch = rep(c("b1", "b2", "b1"), 2),
+        stringsAsFactors = FALSE
+    )
+    rownames(sample_df) <- sample_df$alias
+
+    dds <- DESeq2::DESeqDataSetFromMatrix(
+        countData = SummarizedExperiment::assay(gene_se, "counts"),
+        colData = sample_df,
+        design = ~ batch + condition
+    )
+    out_dir <- tempfile("deseq-fallback-")
+    dir.create(out_dir)
+
+    result <- suppressWarnings(de_run_deseq_with_fallback(
+        dds = dds,
+        contrast_name = "condition_treated_vs_control",
+        out_dir = out_dir
+    ))
+
+    testthat::expect_true(!is.null(result$dds))
+    testthat::expect_true(result$deseq2_dispersion_fallback$method_used %in% c(
+        "parametric",
+        "gene-wise"
+    ))
+    testthat::expect_true(is.logical(result$deseq2_dispersion_fallback$applied))
+    if (isTRUE(result$deseq2_dispersion_fallback$applied)) {
+        testthat::expect_true(file.exists(file.path(
+            out_dir,
+            result$deseq2_dispersion_fallback$diagnostic_file
+        )))
+    }
+})
+
+testthat::test_that("de_run_deseq_with_fallback rethrows non-recoverable errors", {
+    testthat::skip_if_not_installed("DESeq2")
+
+    out_dir <- tempfile("deseq-fallback-error-")
+    dir.create(out_dir)
+
+    testthat::expect_error(
+        de_run_deseq_with_fallback(
+            dds = list(not = "a DESeqDataSet"),
+            contrast_name = "condition_treated_vs_control",
+            out_dir = out_dir
         )
     )
 })
 
+testthat::test_that("de_estimate_dispersions_with_fallback reports method metadata", {
+    testthat::skip_if_not_installed("DESeq2")
+
+    gene_se <- make_test_gene_se()
+    sample_df <- data.frame(
+        alias = colnames(gene_se),
+        condition = rep(c("control", "treated"), each = 3),
+        batch = rep(c("b1", "b2", "b1"), 2),
+        stringsAsFactors = FALSE
+    )
+    rownames(sample_df) <- sample_df$alias
+    dds <- DESeq2::DESeqDataSetFromMatrix(
+        countData = SummarizedExperiment::assay(gene_se, "counts"),
+        colData = sample_df,
+        design = ~ batch + condition
+    )
+    dds <- DESeq2::estimateSizeFactors(dds)
+
+    result <- suppressWarnings(suppressMessages(
+        de_estimate_dispersions_with_fallback(dds, "DESeq2")
+    ))
+
+    testthat::expect_true(result$method_used %in% c(
+        "parametric",
+        "local",
+        "mean",
+        "gene-wise"
+    ))
+    testthat::expect_true(is.logical(result$fallback_applied))
+    testthat::expect_true(!is.null(result$object))
+})
+
+testthat::test_that("de_is_recoverable_dexseq_error recognises expected messages", {
+    testthat::expect_true(de_is_recoverable_dexseq_error(
+        "all gene-wise dispersion estimates are within 2 orders of magnitude"
+    ))
+    testthat::expect_true(de_is_recoverable_dexseq_error("model matrix is not full rank"))
+    testthat::expect_true(de_is_recoverable_dexseq_error("replacement has 1 row, data has 0"))
+    testthat::expect_false(de_is_recoverable_dexseq_error("random unrelated failure"))
+})
+
+testthat::test_that("de_run_dexseq_result records rank-deficiency covariate drops", {
+    testthat::skip_if_not_installed("DESeq2")
+    testthat::skip_if_not_installed("DEXSeq")
+
+    tx_se <- make_test_tx_se(
+        sample_names = c(
+            "control_rep1", "control_rep2", "control_rep3",
+            "treated_rep1", "treated_rep2", "treated_rep3"
+        )
+    )
+    tx_counts <- SummarizedExperiment::assay(tx_se, "counts")
+    tx_meta <- as.data.frame(SummarizedExperiment::rowData(tx_se))
+    coldata <- data.frame(
+        alias = colnames(tx_counts),
+        condition = rep(c("control", "treated"), each = 3),
+        batch = rep(c("control", "treated"), each = 3),
+        stringsAsFactors = FALSE
+    )
+
+    seen_messages <- character(0)
+    result <- withCallingHandlers(
+        tryCatch(
+            de_run_dexseq_result(
+                tx_counts,
+                tx_meta,
+                coldata,
+                condition_column = "condition",
+                covariates = c("batch")
+            ),
+            error = function(err) err
+        ),
+        message = function(m) {
+            seen_messages <<- c(seen_messages, conditionMessage(m))
+            invokeRestart("muffleMessage")
+        }
+    )
+
+    testthat::expect_true(any(grepl(
+        "retrying without it",
+        seen_messages,
+        fixed = TRUE
+    )))
+
+    if (inherits(result, "error")) {
+        testthat::expect_match(
+            conditionMessage(result),
+            "model matrix is not full rank",
+            fixed = TRUE
+        )
+    } else {
+        testthat::expect_equal(result$dexseq_covariates_dropped, "batch")
+        testthat::expect_true(result$dexseq_dispersion_method %in% c(
+            "parametric",
+            "local",
+            "mean",
+            "gene-wise"
+        ))
+        testthat::expect_true(nrow(as.data.frame(result$dxr)) > 0)
+    }
+})
+
 ###
 # Contrast planning and output writing
 #
-# Mock DESeq2/DRIMSeq to avoid slow runtime and test workflow logic:
+# Use small synthetic fixtures with the real DESeq2/DEXSeq path to verify
+# workflow-level planning and output writing:
 # - One contrast created per non-reference condition level (treated vs control, treated2 vs control)
 # - Each contrast subsets to only reference + target samples (not all samples)
 # - Output directories created with correct naming
@@ -409,49 +557,14 @@ testthat::test_that("underspecified designs rejected", {
 # Multi-level design expands to multiple pairwise contrasts (all vs reference).
 # Each contrast subsets samples to just reference + target level.
 testthat::test_that("contrasts expanded and samples subsetted", {
+    testthat::skip_if_not_installed("DESeq2")
+    testthat::skip_if_not_installed("DEXSeq")
+
     fixture_dir <- tempfile("de-multi-")
     dir.create(fixture_dir)
 
     levels <- c("control", "treated", "treated2")
     bundle <- write_de_fixture_bundle(fixture_dir, levels = levels)
-    calls <- new.env(parent = emptyenv())
-    calls$deseq <- list()
-
-    fake_deseq <- function(
-        count_mat,
-        coldata,
-        target_level,
-        reference_level,
-        condition_column,
-        covariates,
-        out_dir,
-        contrast_name
-    ) {
-        calls$deseq[[target_level]] <- coldata$alias
-        result <- data.frame(
-            baseMean = seq_len(nrow(count_mat)),
-            log2FoldChange = rep(1, nrow(count_mat)),
-            lfcSE = rep(0.1, nrow(count_mat)),
-            stat = rep(1, nrow(count_mat)),
-            pvalue = rep(0.05, nrow(count_mat)),
-            padj = rep(0.05, nrow(count_mat)),
-            row.names = rownames(count_mat)
-        )
-        list(dds = structure(list(), class = "fake_dds"), result = result)
-    }
-
-    fake_dexseq <- function(tx_counts, tx_meta, coldata, condition_column, covariates) {
-        dxr <- data.frame(
-            featureID = tx_meta$TXNAME,
-            groupID = tx_meta$GENEID,
-            log2fold = rep(0.5, nrow(tx_meta)),
-            pvalue = rep(0.5, nrow(tx_meta)),
-            padj = rep(0.5, nrow(tx_meta)),
-            exonBaseMean = rep(10, nrow(tx_meta)),
-            row.names = tx_meta$TXNAME
-        )
-        list(dxd = structure(list(), class = "fake_dxd"), dxr = dxr)
-    }
 
     argv <- c(
         bundle,
@@ -463,29 +576,30 @@ testthat::test_that("contrasts expanded and samples subsetted", {
         )
     )
 
-    main_run_de_analysis(
+    suppressWarnings(suppressMessages(main_run_de_analysis(argv)))
-        argv,
-        deseq_runner = fake_deseq,
-        dexseq_runner = fake_dexseq,
-        pdf_fn = function(path) file.create(path),
-        dev_off_fn = function() NULL,
-        plot_ma_fn = function(...) NULL,
-        plot_disp_fn = function(...) NULL,
-        per_gene_q_fn = function(dxr) stats::setNames(c(0.2, 0.3), c("gene1", "gene2")),
-        placeholder_pdf_fn = function(path, label) file.create(path)
-    )
 
     treated_dir <- file.path(argv$out_dir, "condition_treated_vs_control")
     treated2_dir <- file.path(argv$out_dir, "condition_treated2_vs_control")
 
+    treated_samples <- utils::read.delim(
+        file.path(treated_dir, "samples_used.tsv"),
+        check.names = FALSE,
+        stringsAsFactors = FALSE
+    )
+    treated2_samples <- utils::read.delim(
+        file.path(treated2_dir, "samples_used.tsv"),
+        check.names = FALSE,
+        stringsAsFactors = FALSE
+    )
+
     testthat::expect_true(dir.exists(treated_dir))
     testthat::expect_true(dir.exists(treated2_dir))
     testthat::expect_equal(
-        sort(calls$deseq$treated),
+        sort(treated_samples$alias),
         sort(c("control_rep1", "control_rep2", "control_rep3", "treated_rep1", "treated_rep2", "treated_rep3"))
     )
     testthat::expect_equal(
-        sort(calls$deseq$treated2),
+        sort(treated2_samples$alias),
         sort(c("control_rep1", "control_rep2", "control_rep3", "treated2_rep1", "treated2_rep2", "treated2_rep3"))
     )
 })
@@ -493,6 +607,9 @@ testthat::test_that("contrasts expanded and samples subsetted", {
 # Transcript IDs may contain '|' (e.g., Ensembl IDs like ENST00000123.4|ENSG00000456.7).
 # TSV reading defaults to using '|' as separator - verify workflow preserves these IDs.
 testthat::test_that("pipe characters in transcript IDs preserved", {
+    testthat::skip_if_not_installed("DESeq2")
+    testthat::skip_if_not_installed("DEXSeq")
+
     fixture_dir <- tempfile("de-pipe-ids-")
     dir.create(fixture_dir)
     bundle <- write_de_fixture_bundle(fixture_dir, levels = c("control", "treated"))
@@ -508,42 +625,6 @@ testthat::test_that("pipe characters in transcript IDs preserved", {
     S4Vectors::mcols(SummarizedExperiment::rowRanges(tx_se))$TXNAME <- pipe_ids
     saveRDS(tx_se, bundle$transcript_rds)
 
-    captured <- new.env(parent = emptyenv())
-    fake_deseq <- function(
-        count_mat,
-        coldata,
-        target_level,
-        reference_level,
-        condition_column,
-        covariates,
-        out_dir,
-        contrast_name
-    ) {
-        result <- data.frame(
-            baseMean = seq_len(nrow(count_mat)),
-            log2FoldChange = rep(1, nrow(count_mat)),
-            lfcSE = rep(0.1, nrow(count_mat)),
-            stat = rep(1, nrow(count_mat)),
-            pvalue = rep(0.05, nrow(count_mat)),
-            padj = rep(0.05, nrow(count_mat)),
-            row.names = rownames(count_mat)
-        )
-        list(dds = structure(list(), class = "fake_dds"), result = result)
-    }
-    fake_dexseq <- function(tx_counts, tx_meta, coldata, condition_column, covariates) {
-        captured$feature_ids <- tx_meta$TXNAME
-        dxr <- data.frame(
-            featureID = tx_meta$TXNAME,
-            groupID = tx_meta$GENEID,
-            log2fold = rep(0.5, nrow(tx_meta)),
-            pvalue = rep(0.5, nrow(tx_meta)),
-            padj = rep(0.5, nrow(tx_meta)),
-            exonBaseMean = rep(10, nrow(tx_meta)),
-            row.names = tx_meta$TXNAME
-        )
-        list(dxd = structure(list(), class = "fake_dxd"), dxr = dxr)
-    }
-
     argv <- c(
         bundle,
         list(
@@ -554,17 +635,7 @@ testthat::test_that("pipe characters in transcript IDs preserved", {
         )
     )
 
-    main_run_de_analysis(
+    suppressWarnings(suppressMessages(main_run_de_analysis(argv)))
-        argv,
-        deseq_runner = fake_deseq,
-        dexseq_runner = fake_dexseq,
-        pdf_fn = function(path) file.create(path),
-        dev_off_fn = function() NULL,
-        plot_ma_fn = function(...) NULL,
-        plot_disp_fn = function(...) NULL,
-        per_gene_q_fn = function(dxr) stats::setNames(c(0.2, 0.3), c("gene1", "gene2")),
-        placeholder_pdf_fn = function(path, label) file.create(path)
-    )
 
     contrast_dir <- file.path(argv$out_dir, "condition_treated_vs_control")
     dtu_tx <- utils::read.delim(
@@ -577,10 +648,23 @@ testthat::test_that("pipe characters in transcript IDs preserved", {
         check.names = FALSE,
         stringsAsFactors = FALSE
     )
+    de_qc <- jsonlite::read_json(
+        file.path(argv$out_dir, "de_qc_stats.json"),
+        simplifyVector = TRUE
+    )
+    contrast_qc <- de_qc$contrasts[["condition_treated_vs_control"]]
 
-    testthat::expect_equal(captured$feature_ids, pipe_ids)
+    if (identical(contrast_qc$dtu_status, "SUCCESS")) {
         testthat::expect_equal(dtu_tx$featureID, pipe_ids)
         testthat::expect_equal(dexseq$featureID, pipe_ids)
+    } else {
+        testthat::expect_true(file.exists(file.path(
+            contrast_dir,
+            "DTU_ANALYSIS_FAILED.txt"
+        )))
+        testthat::expect_true(nrow(dtu_tx) == 0)
+        testthat::expect_true(nrow(dexseq) == 0)
+    }
 })
 
 ###
@@ -625,9 +709,19 @@ testthat::test_that("CLI integration produces expected outputs", {
     dge <- utils::read.delim(file.path(contrast_dir, "results_dge.tsv"), check.names = FALSE)
     dtu_tx <- utils::read.delim(file.path(contrast_dir, "results_dtu_transcript.tsv"), check.names = FALSE)
     dexseq <- utils::read.delim(file.path(contrast_dir, "results_dexseq.tsv"), check.names = FALSE)
+    de_qc <- jsonlite::read_json(
+        file.path(out_dir, "de_qc_stats.json"),
+        simplifyVector = TRUE
+    )
 
     testthat::expect_gt(nrow(dge), 0)
     testthat::expect_true(all(c("GENEID", "log2FoldChange", "padj") %in% names(dge)))
     testthat::expect_true(all(c("featureID", "groupID", "padj") %in% names(dtu_tx)))
     testthat::expect_true(all(c("featureID", "groupID", "padj") %in% names(dexseq)))
+    testthat::expect_true("analysis_fallbacks" %in% names(de_qc))
+    testthat::expect_true("contrasts" %in% names(de_qc))
+    contrast_qc <- de_qc$contrasts[["condition_treated_vs_control"]]
+    testthat::expect_true("deseq2_dispersion_fallback" %in% names(contrast_qc))
+    testthat::expect_true("dexseq_dispersion_method" %in% names(contrast_qc))
+    testthat::expect_true("dexseq_covariates_dropped" %in% names(contrast_qc))
 })

docs/08_outputs.md

@@ -30,6 +30,12 @@ Output files may be aggregated including information for all samples or provided
 | Differential transcript usage gene summary | de_analysis/{{ contrast }}/results_dtu_gene.tsv | Gene-level DTU summary for one contrast. | aggregated |
 | DEXSeq results | de_analysis/{{ contrast }}/results_dexseq.tsv | Full DEXSeq result table for one contrast. | aggregated |
 | Differential transcript usage plots | de_analysis/{{ contrast }}/results_dtu.pdf | PDF plots generated during DEXSeq analysis for one contrast. | aggregated |
+| Differential analysis QC summary | de_analysis/de_qc_stats.json | Structured DE/DTU QC summary. Use analysis_fallbacks for aggregate counts, and each contrast's deseq2_dispersion_fallback, dexseq_dispersion_method, and dexseq_covariates_dropped fields for interpretation. | aggregated |
+| Differential analysis text summary | de_analysis/de_overall_summary.txt | Human-readable DE/DTU run summary across all contrasts. | aggregated |
+| Per-contrast QC summary | de_analysis/{{ contrast }}/contrast_qc_summary.txt | Human-readable per-contrast DE/DTU QC summary including sample counts and key significance totals. | aggregated |
+| DESeq2 fallback diagnostic | de_analysis/DESeq2_dispersion_fallback_{{ contrast }}.txt | Diagnostic details when DESeq2 falls back to gene-wise dispersion estimation. | aggregated |
+| DTU failure diagnostic | de_analysis/{{ contrast }}/DTU_ANALYSIS_FAILED.txt | Diagnostic details when DEXSeq fails for a contrast. | aggregated |
+| Multiple-testing warning | de_analysis/MULTIPLE_TESTING_WARNING.txt | Family-wise error-rate note generated when multiple contrasts are tested. | aggregated |
 | IGV configuration | igv.json | JSON configuration for viewing the aligned BAMs in IGV. | aggregated |
 | Reference FASTA index | igv_reference/{{ ref_genome_file }}.fai | FAI index for the reference genome published for IGV. | aggregated |
 | Reference GZI index | igv_reference/{{ ref_genome_file }}.gzi | GZI index for a compressed reference genome published for IGV. | aggregated |

output_definition.json

@@ -224,6 +224,54 @@
       "optional": true,
       "type": "aggregated"
     },
+    "de-qc-stats": {
+      "filepath": "de_analysis/de_qc_stats.json",
+      "title": "Differential analysis QC summary",
+      "description": "Structured DE/DTU QC summary. Use analysis_fallbacks for aggregate counts, and each contrast's deseq2_dispersion_fallback, dexseq_dispersion_method, and dexseq_covariates_dropped fields for interpretation.",
+      "mime-type": "application/json",
+      "optional": true,
+      "type": "aggregated"
+    },
+    "de-overall-summary": {
+      "filepath": "de_analysis/de_overall_summary.txt",
+      "title": "Differential analysis text summary",
+      "description": "Human-readable DE/DTU run summary across all contrasts.",
+      "mime-type": "text/plain",
+      "optional": true,
+      "type": "aggregated"
+    },
+    "de-contrast-qc-summary": {
+      "filepath": "de_analysis/{{ contrast }}/contrast_qc_summary.txt",
+      "title": "Per-contrast QC summary",
+      "description": "Human-readable per-contrast DE/DTU QC summary including sample counts and key significance totals.",
+      "mime-type": "text/plain",
+      "optional": true,
+      "type": "aggregated"
+    },
+    "deseq2-dispersion-fallback-diagnostic": {
+      "filepath": "de_analysis/DESeq2_dispersion_fallback_{{ contrast }}.txt",
+      "title": "DESeq2 fallback diagnostic",
+      "description": "Diagnostic details when DESeq2 falls back to gene-wise dispersion estimation.",
+      "mime-type": "text/plain",
+      "optional": true,
+      "type": "aggregated"
+    },
+    "dtu-analysis-failed-diagnostic": {
+      "filepath": "de_analysis/{{ contrast }}/DTU_ANALYSIS_FAILED.txt",
+      "title": "DTU failure diagnostic",
+      "description": "Diagnostic details when DEXSeq fails for a contrast.",
+      "mime-type": "text/plain",
+      "optional": true,
+      "type": "aggregated"
+    },
+    "multiple-testing-warning": {
+      "filepath": "de_analysis/MULTIPLE_TESTING_WARNING.txt",
+      "title": "Multiple-testing warning",
+      "description": "Family-wise error-rate note generated when multiple contrasts are tested.",
+      "mime-type": "text/plain",
+      "optional": true,
+      "type": "aggregated"
+    },
     "igv-config": {
       "filepath": "igv.json",
       "title": "IGV configuration",