diff --git a/README.md b/README.md index 67fcf3e..507b086 100644 --- a/README.md +++ b/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 | diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index d089f6e..0651ab5 100644 --- a/bin/workflow_glue/report.py +++ b/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 - - if ( + sample_size_warnings = _as_string_list( 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 - dispersion_fallbacks = [] - for contrast_name, contrast_data in de_qc.get( - "contrasts", {} - ).items(): - dispersion_file = ( - Path(args.de_dir) - / f"DESeq2_dispersion_fallback_{contrast_name}.txt" - ) - if dispersion_file.exists(): - dispersion_fallbacks.append(contrast_name) - - if dispersion_fallbacks: + if deseq2_gene_wise or dexseq_gene_wise: + gene_wise_details = [] + if deseq2_gene_wise: + gene_wise_details.append( + "DESeq2: " + ", ".join(sorted(deseq2_gene_wise)) + ) + if dexseq_gene_wise: + gene_wise_details.append( + "DEXSeq: " + ", ".join(sorted(dexseq_gene_wise)) + ) _create_warning_banner( - "Dispersion Estimation Fallback: " - f"{len(dispersion_fallbacks)} contrast(s) used " - "gene-wise dispersion (reduced power). " - f"Affected: {', '.join(dispersion_fallbacks)}", + "Gene-wise dispersion fallback used (reduced power). " + + " ".join(gene_wise_details), + level="warning", + ) + 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_value = ( - ", ".join(covariates) - if de_qc.get("covariates") != "none" - else "none" - ) + covariates = _as_string_list(de_qc.get("covariates")) + covariates_value = ", ".join(covariates) if covariates 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 ( - de_qc.get("sample_size_warnings") - and de_qc["sample_size_warnings"] != "none" - ): + if sample_size_warnings: 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" ), } diff --git a/bin/workflow_glue/tests/test_report.py b/bin/workflow_glue/tests/test_report.py index 8c58242..92f2e18 100644 --- a/bin/workflow_glue/tests/test_report.py +++ b/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() diff --git a/bin/workflow_glue_r/R/de_analysis.R b/bin/workflow_glue_r/R/de_analysis.R index 08a2a88..a3b9c97 100644 --- a/bin/workflow_glue_r/R/de_analysis.R +++ b/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) { - tryCatch( +de_extract_disp_gene_est <- function(object) { + 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: %.3f to %.3f", - min(DESeq2::dispersions(dds)), - max(DESeq2::dispersions(dds)) - ), + dispersion_range, "", "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 + object <- de_set_dispersions(object, de_extract_disp_gene_est(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( - 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 -) { +main_run_de_analysis <- function(argv) { 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")) - plot_ma_fn(dge_run$result) - dev_off_fn() + grDevices::pdf(file.path(contrast_dir, "results_dge.pdf")) + DESeq2::plotMA(dge_run$result) + 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")) - plot_ma_fn(dex_res$dxr, cex = 0.8, alpha = 0.05) - plot_disp_fn(dex_res$dxd) - dev_off_fn() + grDevices::pdf(file.path(contrast_dir, "results_dtu.pdf")) + DESeq2::plotMA(dex_res$dxr, cex = 0.8, alpha = 0.05) + DESeq2::plotDispEsts(dex_res$dxd) + 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"), diff --git a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R b/bin/workflow_glue_r/tests/testthat/test_de_analysis.R index eb31a46..de13f5b 100644 --- a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R +++ b/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( - argv, - deseq_runner = function(...) stop("runner should not be called"), - dexseq_runner = function(...) stop("runner should not be called") - ), + main_run_de_analysis(argv), "Transcript rowData must contain GENEID" ) }) @@ -347,11 +343,7 @@ testthat::test_that("underspecified designs rejected", { ) testthat::expect_error( - suppressWarnings(main_run_de_analysis( - argv, - deseq_runner = function(...) stop("runner should not be called"), - dexseq_runner = function(...) stop("runner should not be called") - )), + suppressWarnings(main_run_de_analysis(argv)), "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( - argv, - deseq_runner = function(...) stop("runner should not be called"), - dexseq_runner = function(...) stop("runner should not be called") - ), + main_run_de_analysis(argv), "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( - argv, - deseq_runner = function(...) stop("runner should not be called"), - dexseq_runner = function(...) stop("runner should not be called") + main_run_de_analysis(argv) + ) +}) + +### +# 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( - 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) - ) + suppressWarnings(suppressMessages(main_run_de_analysis(argv))) 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( - 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) - ) + suppressWarnings(suppressMessages(main_run_de_analysis(argv))) 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) - testthat::expect_equal(dtu_tx$featureID, pipe_ids) - testthat::expect_equal(dexseq$featureID, 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)) }) diff --git a/docs/08_outputs.md b/docs/08_outputs.md index bf9de17..3729e61 100644 --- a/docs/08_outputs.md +++ b/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 | diff --git a/output_definition.json b/output_definition.json index 88895e4..f57a98e 100644 --- a/output_definition.json +++ b/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",