diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 929a88f..82eaf6c 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -65,7 +65,7 @@ docker-run: parallel: matrix: - MATRIX_NAME: [ - "int_discover_dna", "int_fixed_rna", + "int_discover_dna", "int_fixed_rna", "int_de_control_vs_control", "smoke_discover", "smoke_fixed", "smoke_direct_rna", "smoke_de", ] rules: @@ -77,7 +77,7 @@ docker-run: when: never # Integration: larger discover-mode run on representative cDNA test bundle. - - if: $MATRIX_NAME == "discover" + - if: $MATRIX_NAME == "int_discover_dna" variables: NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;" NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq \ @@ -102,6 +102,37 @@ docker-run: NF_IGNORE_PROCESSES: > gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,decompress_ref, build_minimap_index,validate_ref_annotation,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome + + # Integration: demo-data control-vs-control to test pathological DE/DTU failure modes and graceful failure handling. + # See also unit tests in test_de_analysis.R + - if: $MATRIX_NAME == "int_de_control_vs_control" + variables: + NF_BEFORE_SCRIPT: "\ + mkdir -p ${CI_PROJECT_NAME}/data/ \ + && echo 'Downlading demo data bundle and config' \ + && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz \ + && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ \ + && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config \ + && echo 'Simulating control-vs-control by copying each sample to a replicate with identical reads' \ + && rm -rf ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode04 ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode05 ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode06 \ + && cp -R ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode01 ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode04 \ + && cp -R ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode02 ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode05 \ + && cp -R ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode03 ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq/barcode06 \ + && echo 'Finished data preparation step' \ + " + NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \ + --de_analysis \ + --ref_genome ${CI_PROJECT_NAME}/data/differential_expression/hg38_chr20.fa \ + --ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gtf \ + --direct_rna \ + --sample_sheet test_data/sample_sheet.csv \ + -c ${CI_PROJECT_NAME}/data/demo.nextflow.config" + AFTER_NEXTFLOW_CMD: > + test -f ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json && + test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dge.tsv && + test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dtu_transcript.tsv && + jq -e '.contrasts["condition_treated_vs_control"].dge_status | IN("SUCCESS", "FAILED")' ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json >/dev/null && + jq -e '.contrasts["condition_treated_vs_control"].dtu_status | IN("SUCCESS", "FAILED")' ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json >/dev/null # Smoke: quick discover-mode sanity check for core cohort and per-sample outputs. - if: $MATRIX_NAME == "smoke_discover" diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index f13f1cb..af40af4 100644 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -177,6 +177,7 @@ def _collect_de_method_rows(de_qc): deseq2_gene_wise = [] dexseq_gene_wise = [] dexseq_covariate_drops = [] + failed_dge = [] for contrast_name, contrast_data in de_qc.get("contrasts", {}).items(): fallback = contrast_data.get("deseq2_dispersion_fallback") or {} @@ -201,6 +202,9 @@ def _collect_de_method_rows(de_qc): if dropped_covariates: dexseq_covariate_drops.append((contrast_name, dropped_covariates)) + if contrast_data.get("dge_status") == "FAILED": + failed_dge.append(contrast_name) + rows.append( { "Contrast": contrast_name, @@ -215,11 +219,12 @@ def _collect_de_method_rows(de_qc): "DEXSeq covariates dropped": ( ", ".join(dropped_covariates) if dropped_covariates else "none" ), + "DGE status": contrast_data.get("dge_status", "N/A"), "DTU status": contrast_data.get("dtu_status", "N/A"), } ) - return rows, deseq2_gene_wise, dexseq_gene_wise, dexseq_covariate_drops + return rows, deseq2_gene_wise, dexseq_gene_wise, dexseq_covariate_drops, failed_dge def main(args): @@ -529,6 +534,7 @@ def main(args): deseq2_gene_wise, dexseq_gene_wise, dexseq_covariate_drops, + failed_dge, ) = _collect_de_method_rows(de_qc) if sample_size_warnings: @@ -578,6 +584,17 @@ def main(args): ) has_warnings = True + if failed_dge: + _create_warning_banner( + "DGE Analysis Failed: " + f"{len(failed_dge)} contrast(s) could not complete " + "DGE testing. " + f"Affected: {', '.join(failed_dge)}. " + "See DGE_ANALYSIS_FAILED.txt files for details.", + level="danger", + ) + has_warnings = True + # Check for failed DTU analyses failed_dtu = [] for contrast_name, contrast_data in de_qc.get( @@ -681,16 +698,25 @@ def main(args): f"vs " f"{contrast_data.get('n_reference', 0)}" ), + "DGE Status": contrast_data.get( + "dge_status", "N/A" + ), "DGE Significant (FDR<0.05)": ( contrast_data.get( "dge_significant_fdr05", 0 ) + if contrast_data.get("dge_status") == "SUCCESS" + else "N/A" ), - "DGE Up": contrast_data.get( - "dge_upregulated", 0 + "DGE Up": ( + contrast_data.get("dge_upregulated", 0) + if contrast_data.get("dge_status") == "SUCCESS" + else "N/A" ), - "DGE Down": contrast_data.get( - "dge_downregulated", 0 + "DGE Down": ( + contrast_data.get("dge_downregulated", 0) + if contrast_data.get("dge_status") == "SUCCESS" + else "N/A" ), "DTU Status": contrast_data.get( "dtu_status", "N/A" @@ -751,6 +777,15 @@ def main(args): ), } ) + if failed_dge: + warnings_data.append( + { + "Warning Type": "DGE Failure", + "Details": ( + f"{len(failed_dge)} contrasts failed" + ), + } + ) warnings_df = pd.DataFrame(warnings_data) DataTable.from_pandas( @@ -766,7 +801,19 @@ def main(args): # Check for contrast-specific warnings if de_qc and contrast in de_qc.get("contrasts", {}): contrast_data = de_qc["contrasts"][contrast] - if contrast_data.get("dtu_power_warning"): + if contrast_data.get("dge_status") == "FAILED": + _create_warning_banner( + "DGE analysis failed for this contrast. " + f"See {contrast}/" + "DGE_ANALYSIS_FAILED.txt for detailed " + "explanation.", + level="danger", + ) + p( + "Empty results indicate analysis failure, " + "not 'no DGE detected'." + ) + elif contrast_data.get("dtu_power_warning"): with div( style=( "padding: 10px; margin-bottom: 10px; " @@ -801,6 +848,9 @@ def main(args): "explanation.", level="danger", ) + failure_hint = contrast_data.get("dtu_failure_hint") + if failure_hint: + p(f"Probable cause: {failure_hint}") p( "Empty results indicate analysis failure, " "not 'no DTU detected'." diff --git a/bin/workflow_glue_r/R/de_analysis.R b/bin/workflow_glue_r/R/de_analysis.R index 0205f92..7b10b22 100644 --- a/bin/workflow_glue_r/R/de_analysis.R +++ b/bin/workflow_glue_r/R/de_analysis.R @@ -201,29 +201,64 @@ de_run_deseq_with_fallback <- function( ) de_out <- tryCatch( - DESeq2::DESeq(dds, quiet = TRUE, sfType = sf_type), + list( + dds = DESeq2::DESeq(dds, quiet = TRUE, sfType = sf_type), + deseq2_dispersion_fallback = fallback_info + ), error = function(err) { - if (!grepl( - "all gene-wise dispersion estimates are within 2 orders of magnitude", - conditionMessage(err), + msg <- conditionMessage(err) + zero_geometric_means <- grepl( + "every gene contains at least one zero", + msg, fixed = TRUE - )) { + ) + + # Size factor estimation failure: all genes have at least one zero. + # de_choose_size_factor_type should prevent this proactively, but as + # a safety net switch to poscounts so the gene-wise fallback below + # can call estimateSizeFactors without failing again. + fallback_sf_type <- if (zero_geometric_means) "poscounts" else sf_type + + recoverable <- grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + msg, + fixed = TRUE + ) || grepl( + "newsplit: out of vertex space", + msg, + fixed = TRUE + ) || grepl( + "every gene contains at least one zero", + msg, + fixed = TRUE + ) + if (!recoverable) { stop(err) } - warning( - "STATISTICAL POWER REDUCED: DESeq2 dispersion estimation failed for ", - contrast_name, - ".\n", - "This usually indicates:\n", - " 1. Too few replicates (recommend n>=3 per group)\n", - " 2. High biological variability\n", - " 3. Poor data quality\n", - "Falling back to gene-wise dispersion (no information sharing).\n", - "Results will have reduced power and wider confidence intervals." - ) + if (zero_geometric_means) { + warning( + "STATISTICAL POWER REDUCED: DESeq2 default size-factor estimation failed for ", + contrast_name, + " because every gene contains at least one zero.\n", + "Switching to sfType='poscounts' and using gene-wise dispersion estimates.\n", + "Results will have reduced power and wider confidence intervals." + ) + } else { + warning( + "STATISTICAL POWER REDUCED: DESeq2 dispersion estimation failed for ", + contrast_name, + ".\n", + "This usually indicates:\n", + " 1. Too few replicates (recommend n>=3 per group)\n", + " 2. High biological variability\n", + " 3. Poor data quality\n", + "Falling back to gene-wise dispersion (no information sharing).\n", + "Results will have reduced power and wider confidence intervals." + ) + } - dds <- DESeq2::estimateSizeFactors(dds, type = sf_type) + dds <- DESeq2::estimateSizeFactors(dds, type = fallback_sf_type) dds <- DESeq2::estimateDispersionsGeneEst(dds) dispersions_setter <- get("dispersions<-", envir = asNamespace("DESeq2")) dds <- dispersions_setter(dds, value = S4Vectors::mcols(dds)$dispGeneEst) @@ -241,17 +276,35 @@ de_run_deseq_with_fallback <- function( } diag_content <- c( - "DESeq2 Dispersion Estimation Fallback Applied", - "==============================================", + if (zero_geometric_means) { + "DESeq2 Poscounts and Gene-wise Fallback Applied" + } else { + "DESeq2 Dispersion Estimation Fallback Applied" + }, + if (zero_geometric_means) { + "=============================================" + } else { + "==============================================" + }, "", sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), sprintf("Contrast: %s", contrast_name), sprintf("Samples: %d", ncol(dds)), sprintf("Genes tested: %d", nrow(dds)), + sprintf("Size factor method: %s", fallback_sf_type), dispersion_range, "", "WHAT HAPPENED:", - " Curve fitting failed. Using gene-wise dispersion estimates.", + if (zero_geometric_means) { + " Every gene contained at least one zero, so DESeq2's default size-factor estimator was undefined." + } else { + " Curve fitting failed. Using gene-wise dispersion estimates." + }, + if (zero_geometric_means) { + " The analysis switched to sfType='poscounts' and used gene-wise dispersion estimates." + } else { + NULL + }, "", "IMPLICATIONS:", " - No information sharing across genes", @@ -260,14 +313,34 @@ de_run_deseq_with_fallback <- function( " - More conservative results (fewer discoveries)", "", "LIKELY CAUSES:", - " 1. Too few replicates (recommend n>=3 per group)", - " 2. High biological variability", - " 3. Poor data quality or outlier samples", + if (zero_geometric_means) { + c( + " 1. Sparse count matrices with many zeros", + " 2. Low-expression features dominating the design", + " 3. Small or weakly separated groups" + ) + } else { + c( + " 1. Too few replicates (recommend n>=3 per group)", + " 2. High biological variability", + " 3. Poor data quality or outlier samples" + ) + }, "", "RECOMMENDATIONS:", - " - Add more biological replicates if possible", - " - Check sample quality metrics", - " - Consider filtering low-count genes more stringently" + if (zero_geometric_means) { + c( + " - Check whether the contrast is dominated by zeros or nearly identical samples", + " - Add more biological replicates if possible", + " - Consider filtering very low-count genes more stringently" + ) + } else { + c( + " - Add more biological replicates if possible", + " - Check sample quality metrics", + " - Consider filtering low-count genes more stringently" + ) + } ) diag_file <- file.path( @@ -278,18 +351,19 @@ de_run_deseq_with_fallback <- function( ) ) writeLines(diag_content, diag_file) - fallback_info <<- list( - applied = TRUE, - method_used = "gene-wise", - reason = conditionMessage(err), - diagnostic_file = basename(diag_file), - size_factor_type = sf_type + list( + dds = DESeq2::nbinomWaldTest(dds), + deseq2_dispersion_fallback = list( + applied = TRUE, + method_used = "gene-wise", + reason = conditionMessage(err), + diagnostic_file = basename(diag_file), + size_factor_type = fallback_sf_type + ) ) - - DESeq2::nbinomWaldTest(dds) } ) - list(dds = de_out, deseq2_dispersion_fallback = fallback_info) + de_out } de_estimate_dispersions_with_fallback <- function( @@ -297,6 +371,15 @@ de_estimate_dispersions_with_fallback <- function( context_label, allow_gene_est = TRUE ) { + de_dispersion_recoverable <- function(msg) { + grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + msg, fixed = TRUE + ) || grepl( + "newsplit: out of vertex space", + msg, fixed = TRUE + ) + } tryCatch( list( object = DESeq2::estimateDispersions(object), @@ -304,11 +387,7 @@ de_estimate_dispersions_with_fallback <- function( fallback_applied = FALSE ), error = function(err) { - if (!grepl( - "all gene-wise dispersion estimates are within 2 orders of magnitude", - conditionMessage(err), - fixed = TRUE - )) { + if (!de_dispersion_recoverable(conditionMessage(err))) { stop(err) } @@ -323,11 +402,7 @@ de_estimate_dispersions_with_fallback <- function( fallback_applied = TRUE ), error = function(local_err) { - if (!grepl( - "all gene-wise dispersion estimates are within 2 orders of magnitude", - conditionMessage(local_err), - fixed = TRUE - )) { + if (!de_dispersion_recoverable(conditionMessage(local_err))) { stop(local_err) } @@ -342,11 +417,7 @@ de_estimate_dispersions_with_fallback <- function( fallback_applied = TRUE ), error = function(mean_err) { - if (!grepl( - "all gene-wise dispersion estimates are within 2 orders of magnitude", - conditionMessage(mean_err), - fixed = TRUE - )) { + if (!de_dispersion_recoverable(conditionMessage(mean_err))) { stop(mean_err) } if (!allow_gene_est) { @@ -389,9 +460,41 @@ de_is_recoverable_dexseq_error <- function(message_text) { "replacement has 1 row, data has 0", message_text, fixed = TRUE + ) || grepl( + "newsplit: out of vertex space", + message_text, + fixed = TRUE + ) || grepl( + "every gene contains at least one zero", + message_text, + fixed = TRUE ) } +de_dexseq_failure_hint <- function(message_text) { + if (grepl("nb of cols in 'assay'", message_text, fixed = TRUE)) { + return( + paste( + "Probable cause: degenerate or near-identical sample groups can", + "drive DEXSeq into an internal object-shape mismatch during model", + "fitting. This is often seen when the contrast has little real", + "separation or too few informative non-zero features." + ) + ) + } + if (grepl("'x' must be an array of at least two dimensions", message_text, fixed = TRUE)) { + return( + paste( + "Probable cause: degenerate or near-identical sample groups can", + "leave DEXSeq with too little informative structure for downstream", + "gene-level DTU summarisation. This is often seen when the contrast", + "has little real separation or too few informative non-zero features." + ) + ) + } + NULL +} + de_write_placeholder_pdf <- function(path, label) { grDevices::pdf(path) graphics::plot.new() @@ -399,6 +502,50 @@ de_write_placeholder_pdf <- function(path, label) { grDevices::dev.off() } +de_write_dtu_failure_outputs <- function( + contrast_dir, + target_level, + reference_level, + contrast_samples, + condition_column, + tx_counts, + message_text, + failure_hint +) { + failure_content <- c( + "DTU Analysis Failed", + "===================", + "", + sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), + sprintf("Contrast: %s vs %s", target_level, reference_level), + sprintf( + "Samples: %d (%d %s, %d %s)", + nrow(contrast_samples), + sum(contrast_samples[[condition_column]] == target_level), + target_level, + sum(contrast_samples[[condition_column]] == reference_level), + reference_level + ), + sprintf("Transcripts: %d", nrow(tx_counts)), + "", + "ERROR MESSAGE:", + sprintf(" %s", message_text), + if (!is.null(failure_hint)) { + c( + "", + "PROBABLE CAUSE:", + sprintf(" %s", failure_hint) + ) + } else { + NULL + }, + "", + "DTU RESULTS CANNOT BE INTERPRETED", + "" + ) + writeLines(failure_content, file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt")) +} + de_run_deseq2_result <- function( count_mat, coldata, @@ -453,9 +600,7 @@ de_run_dexseq_result <- function( for (covariate in covariates) { coldata[[covariate]] <- factor(coldata[[covariate]]) } - dropped_covariates <- character(0) - - run_inner <- function(active_covariates) { + run_inner <- function(active_covariates, dropped_covariates = character(0)) { covariate_exon_terms <- if (length(active_covariates) > 0) { paste0(active_covariates, ":exon") } else { @@ -498,7 +643,8 @@ de_run_dexseq_result <- function( dxd = dxd, dxr = dxr, dexseq_dispersion_method = dispersion_method, - dexseq_size_factor_type = dexseq_sf_type + dexseq_size_factor_type = dexseq_sf_type, + dexseq_covariates_dropped = dropped_covariates ) }, error = function(err) { if (length(active_covariates) == 0 || !grepl( @@ -511,21 +657,23 @@ 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, "'; retrying without it." ) - run_inner(kept_covariates) + run_inner( + kept_covariates, + c(dropped_covariates, dropped_covariate) + ) }) } - result <- run_inner(covariates) - result$dexseq_covariates_dropped <- dropped_covariates - result + run_inner(covariates) } +de_dge_columns <- c("GENEID", "baseMean", "log2FoldChange", "lfcSE", "stat", "pvalue", "padj") + de_dtu_transcript_columns <- c( "featureID", "groupID", @@ -557,6 +705,81 @@ de_extract_dtu_transcript_table <- function(dex_df, contrast_name) { workflow_glue_r_normalise_tsv_df(tx_dtu) } +de_postprocess_dexseq_result <- function( + dex_res, + contrast_dir, + target_level, + reference_level, + contrast_samples, + condition_column, + tx_counts, + contrast_name, + per_gene_qvalue_fn = DEXSeq::perGeneQValue, + plot_writer = function(dex_res, contrast_dir) { + 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() + } +) { + tryCatch( + { + dex_df <- as.data.frame(dex_res$dxr) + dex_df <- workflow_glue_r_normalise_tsv_df(dex_df) + tx_dtu <- de_extract_dtu_transcript_table(dex_df, contrast_name) + gene_q <- per_gene_qvalue_fn(dex_res$dxr) + gene_dtu <- data.frame( + GENEID = names(gene_q), + qval = unname(gene_q), + row.names = NULL + ) + plot_writer(dex_res, contrast_dir) + list( + ok = TRUE, + dex_df = dex_df, + tx_dtu = tx_dtu, + gene_dtu = gene_dtu, + failure_hint = NULL + ) + }, + error = function(err) { + message_text <- conditionMessage(err) + failure_hint <- de_dexseq_failure_hint(message_text) + + warning( + "DEXSeq post-processing failed for contrast ", + target_level, + " vs ", + reference_level, + if (!is.null(failure_hint)) { + paste0("\nProbable cause: ", failure_hint) + } else { + "" + }, + "\nError: ", + message_text + ) + de_write_dtu_failure_outputs( + contrast_dir = contrast_dir, + target_level = target_level, + reference_level = reference_level, + contrast_samples = contrast_samples, + condition_column = condition_column, + tx_counts = tx_counts, + message_text = message_text, + failure_hint = failure_hint + ) + list( + ok = FALSE, + dex_df = NULL, + tx_dtu = NULL, + gene_dtu = NULL, + failure_hint = failure_hint + ) + } + ) +} + main_run_de_analysis <- function(args) { set.seed(42) dir.create(args$out_dir, showWarnings = FALSE, recursive = TRUE) @@ -693,7 +916,9 @@ main_run_de_analysis <- function(args) { ), dexseq_size_factor_method = "ratio", dexseq_dispersion_method = "parametric", - dexseq_covariates_dropped = list() + dexseq_covariates_dropped = list(), + dtu_failure_hint = NULL, + dge_status = "SUCCESS" ) if (nrow(contrast_samples) < 6) { @@ -709,87 +934,25 @@ main_run_de_analysis <- function(args) { contrast_qc$genes_tested <- nrow(gene_counts) contrast_qc$transcripts_tested <- nrow(tx_counts) - dge_run <- de_run_deseq2_result( - gene_counts, - contrast_samples, - target_level, - reference_level, - args$condition_column, - covariates, - args$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" - } - if (!is.null(fallback$size_factor_type) && !identical(fallback$size_factor_type, "")) { - contrast_qc$deseq2_size_factor_method <- fallback$size_factor_type - } - 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) - dge_res <- workflow_glue_r_normalise_tsv_df(dge_res) - - contrast_qc$dge_total_genes <- nrow(dge_res) - contrast_qc$dge_significant_fdr05 <- sum(dge_res$padj < 0.05, na.rm = TRUE) - contrast_qc$dge_significant_fdr01 <- sum(dge_res$padj < 0.01, na.rm = TRUE) - contrast_qc$dge_upregulated <- sum( - dge_res$padj < 0.05 & dge_res$log2FoldChange > 0, - na.rm = TRUE - ) - contrast_qc$dge_downregulated <- sum( - dge_res$padj < 0.05 & dge_res$log2FoldChange < 0, - na.rm = TRUE - ) - - utils::write.table( - dge_res[order(dge_res$padj), ], - file = file.path(contrast_dir, "results_dge.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE - ) - - grDevices::pdf(file.path(contrast_dir, "results_dge.pdf")) - DESeq2::plotMA(dge_run$result) - grDevices::dev.off() - - dex_res <- tryCatch( - de_run_dexseq_result( - tx_counts, - tx_meta, + dge_run <- tryCatch( + de_run_deseq2_result( + gene_counts, contrast_samples, + target_level, + reference_level, args$condition_column, - covariates + covariates, + args$out_dir, + contrast_name ), error = function(err) { - message_text <- conditionMessage(err) - if (!de_is_recoverable_dexseq_error(message_text)) { - stop(err) - } - warning( - "DEXSeq failed for contrast ", - target_level, - " vs ", - reference_level, - "\nError: ", - message_text + "DESeq2 DGE failed for contrast ", + target_level, " vs ", reference_level, + "\nError: ", conditionMessage(err) ) - failure_content <- c( - "DTU Analysis Failed", + "DGE Analysis Failed", "===================", "", sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), @@ -802,19 +965,135 @@ main_run_de_analysis <- function(args) { sum(contrast_samples[[args$condition_column]] == reference_level), reference_level ), - sprintf("Transcripts: %d", nrow(tx_counts)), + sprintf("Genes: %d", nrow(gene_counts)), "", "ERROR MESSAGE:", - sprintf(" %s", message_text), + sprintf(" %s", conditionMessage(err)), "", - "DTU RESULTS CANNOT BE INTERPRETED", + "DGE RESULTS CANNOT BE INTERPRETED", "" ) - writeLines(failure_content, file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt")) + writeLines(failure_content, file.path(contrast_dir, "DGE_ANALYSIS_FAILED.txt")) NULL } ) + if (is.null(dge_run)) { + dge_res <- workflow_glue_r_empty_tsv(de_dge_columns) + contrast_qc$dge_status <- "FAILED" + contrast_qc$dge_total_genes <- 0 + contrast_qc$dge_significant_fdr05 <- 0 + contrast_qc$dge_significant_fdr01 <- 0 + contrast_qc$dge_upregulated <- 0 + contrast_qc$dge_downregulated <- 0 + de_write_placeholder_pdf( + file.path(contrast_dir, "results_dge.pdf"), + "DESeq2 did not converge for this contrast.\nSee DGE_ANALYSIS_FAILED.txt for details." + ) + } else { + 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" + } + if (!is.null(fallback$size_factor_type) && !identical(fallback$size_factor_type, "")) { + contrast_qc$deseq2_size_factor_method <- fallback$size_factor_type + } + 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) + dge_res <- workflow_glue_r_normalise_tsv_df(dge_res) + contrast_qc$dge_status <- "SUCCESS" + contrast_qc$dge_total_genes <- nrow(dge_res) + contrast_qc$dge_significant_fdr05 <- sum(dge_res$padj < 0.05, na.rm = TRUE) + contrast_qc$dge_significant_fdr01 <- sum(dge_res$padj < 0.01, na.rm = TRUE) + contrast_qc$dge_upregulated <- sum( + dge_res$padj < 0.05 & dge_res$log2FoldChange > 0, + na.rm = TRUE + ) + contrast_qc$dge_downregulated <- sum( + dge_res$padj < 0.05 & dge_res$log2FoldChange < 0, + na.rm = TRUE + ) + grDevices::pdf(file.path(contrast_dir, "results_dge.pdf")) + DESeq2::plotMA(dge_run$result) + grDevices::dev.off() + } + + utils::write.table( + dge_res[order(dge_res$padj), ], + file = file.path(contrast_dir, "results_dge.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + + dex_run <- tryCatch( + list( + ok = TRUE, + result = de_run_dexseq_result( + tx_counts, + tx_meta, + contrast_samples, + args$condition_column, + covariates + ), + failure_hint = NULL + ), + error = function(err) { + message_text <- conditionMessage(err) + is_known_failure <- de_is_recoverable_dexseq_error(message_text) + failure_hint <- de_dexseq_failure_hint(message_text) + + warning( + "DEXSeq failed for contrast ", + target_level, + " vs ", + reference_level, + if (!is_known_failure) { + "\nThis error did not match a known statistical fallback pattern." + } else { + "" + }, + if (!is.null(failure_hint)) { + paste0("\nProbable cause: ", failure_hint) + } else { + "" + }, + "\nError: ", + message_text + ) + de_write_dtu_failure_outputs( + contrast_dir = contrast_dir, + target_level = target_level, + reference_level = reference_level, + contrast_samples = contrast_samples, + condition_column = args$condition_column, + tx_counts = tx_counts, + message_text = message_text, + failure_hint = failure_hint + ) + list(ok = FALSE, result = NULL, failure_hint = failure_hint) + } + ) + + if (!isTRUE(dex_run$ok)) { + dex_res <- NULL + dtu_failure_hint <- dex_run$failure_hint + } else { + dex_res <- dex_run$result + dtu_failure_hint <- NULL + } + if (is.null(dex_res)) { dex_df <- workflow_glue_r_empty_tsv(de_dtu_transcript_columns) tx_dtu <- dex_df @@ -824,6 +1103,7 @@ main_run_de_analysis <- function(args) { "DEXSeq did not converge for this contrast.\nSee DTU_ANALYSIS_FAILED.txt for details." ) contrast_qc$dtu_status <- "FAILED" + contrast_qc$dtu_failure_hint <- dtu_failure_hint contrast_qc$dtu_significant_transcripts <- 0 contrast_qc$dtu_significant_genes <- 0 } else { @@ -836,28 +1116,39 @@ main_run_de_analysis <- function(args) { 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 <- de_extract_dtu_transcript_table( - dex_df, - paste(target_level, reference_level, sep = "_") + + postprocess_run <- de_postprocess_dexseq_result( + dex_res = dex_res, + contrast_dir = contrast_dir, + target_level = target_level, + reference_level = reference_level, + contrast_samples = contrast_samples, + condition_column = args$condition_column, + tx_counts = tx_counts, + contrast_name = paste(target_level, reference_level, sep = "_") ) - gene_q <- DEXSeq::perGeneQValue(dex_res$dxr) - gene_dtu <- data.frame( - GENEID = names(gene_q), - qval = unname(gene_q), - row.names = NULL - ) + if (!isTRUE(postprocess_run$ok)) { + dex_df <- workflow_glue_r_empty_tsv(de_dtu_transcript_columns) + tx_dtu <- dex_df + gene_dtu <- workflow_glue_r_empty_tsv(c("GENEID", "qval")) + 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." + ) + contrast_qc$dtu_status <- "FAILED" + contrast_qc$dtu_failure_hint <- postprocess_run$failure_hint + contrast_qc$dtu_significant_transcripts <- 0 + contrast_qc$dtu_significant_genes <- 0 + } else { + dex_df <- postprocess_run$dex_df + tx_dtu <- postprocess_run$tx_dtu + gene_dtu <- postprocess_run$gene_dtu - contrast_qc$dtu_status <- "SUCCESS" - 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) - - 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() + contrast_qc$dtu_status <- "SUCCESS" + 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) + } } utils::write.table( @@ -899,12 +1190,25 @@ main_run_de_analysis <- function(args) { sprintf(" Total samples: %d", contrast_qc$n_samples), "", "DGE Results:", + sprintf(" Status: %s", contrast_qc$dge_status), sprintf(" Genes tested: %d", contrast_qc$genes_tested), sprintf(" Size factor method: %s", contrast_qc$deseq2_size_factor_method), - sprintf(" Significant (FDR < 0.05): %d", contrast_qc$dge_significant_fdr05), - sprintf(" Significant (FDR < 0.01): %d", contrast_qc$dge_significant_fdr01), - sprintf(" Upregulated: %d", contrast_qc$dge_upregulated), - sprintf(" Downregulated: %d", contrast_qc$dge_downregulated), + if (contrast_qc$dge_status == "SUCCESS") { + c( + sprintf(" Significant (FDR < 0.05): %d", contrast_qc$dge_significant_fdr05), + sprintf(" Significant (FDR < 0.01): %d", contrast_qc$dge_significant_fdr01), + sprintf(" Upregulated: %d", contrast_qc$dge_upregulated), + sprintf(" Downregulated: %d", contrast_qc$dge_downregulated) + ) + } else { + c( + " Significant (FDR < 0.05): N/A", + " Significant (FDR < 0.01): N/A", + " Upregulated: N/A", + " Downregulated: N/A", + " See DGE_ANALYSIS_FAILED.txt for details" + ) + }, "", "DTU Results:", sprintf(" Status: %s", contrast_qc$dtu_status), @@ -919,7 +1223,14 @@ main_run_de_analysis <- function(args) { sprintf(" Genes with DTU (q < 0.05): %d", contrast_qc$dtu_significant_genes) ) } else { - " See DTU_ANALYSIS_FAILED.txt for details" + c( + if (!is.null(contrast_qc$dtu_failure_hint)) { + paste0(" Probable cause: ", contrast_qc$dtu_failure_hint) + } else { + NULL + }, + " See DTU_ANALYSIS_FAILED.txt for details" + ) }, if (!is.null(contrast_qc$dtu_power_warning)) paste0(" WARNING: ", contrast_qc$dtu_power_warning) else NULL, "" @@ -928,6 +1239,15 @@ main_run_de_analysis <- function(args) { de_qc_stats$contrasts[[contrast_name]] <- contrast_qc } + failed_dge_contrasts <- names(Filter( + function(cqc) identical(cqc$dge_status, "FAILED"), + de_qc_stats$contrasts + )) + failed_dtu_contrasts <- names(Filter( + function(cqc) identical(cqc$dtu_status, "FAILED"), + de_qc_stats$contrasts + )) + deseq2_dispersion_fallbacks <- names(Filter( function(cqc) isTRUE(cqc$deseq2_dispersion_fallback$applied), de_qc_stats$contrasts @@ -965,6 +1285,10 @@ main_run_de_analysis <- function(args) { integer(1) )) de_qc_stats$analysis_fallbacks <- list( + failed_dge_contrasts = length(failed_dge_contrasts), + failed_dge_contrast_names = as.list(failed_dge_contrasts), + failed_dtu_contrasts = length(failed_dtu_contrasts), + failed_dtu_contrast_names = as.list(failed_dtu_contrasts), 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), @@ -1015,8 +1339,14 @@ main_run_de_analysis <- function(args) { "", sprintf(" %s:", cname), sprintf(" Samples: %d (%d vs %d)", cqc$n_samples, cqc$n_target, cqc$n_reference), - sprintf(" DGE significant: %d genes (FDR<0.05)", cqc$dge_significant_fdr05), + sprintf(" DGE status: %s", cqc$dge_status), + if (cqc$dge_status == "SUCCESS") { + sprintf(" DGE significant: %d genes (FDR<0.05)", cqc$dge_significant_fdr05) + } else { + " DGE significant: N/A (analysis failed)" + }, sprintf(" DESeq2 size factors: %s", cqc$deseq2_size_factor_method), + if (cqc$dge_status != "SUCCESS") " See DGE_ANALYSIS_FAILED.txt for details" else NULL, sprintf(" DTU status: %s", cqc$dtu_status), sprintf(" DEXSeq size factors: %s", cqc$dexseq_size_factor_method), if (cqc$dtu_status == "SUCCESS") sprintf(" DTU significant: %d genes", cqc$dtu_significant_genes) else NULL diff --git a/bin/workflow_glue_r/tests/testthat/helper-fixtures.R b/bin/workflow_glue_r/tests/testthat/helper-fixtures.R index eedfb08..0a7e59f 100644 --- a/bin/workflow_glue_r/tests/testthat/helper-fixtures.R +++ b/bin/workflow_glue_r/tests/testthat/helper-fixtures.R @@ -308,3 +308,86 @@ expect_bam_fixture_built <- function(reference, reads, out_dir, alias = "sampleA bam_path } + +#' Create a degenerate DE fixture where both condition groups have identical +#' zero-heavy counts, simulating the control-vs-control case from the bug reports. +#' +#' Properties: +#' - Every transcript has at least one zero (forces DEXSeq poscounts path) +#' - Every gene has at least one zero (forces DESeq2 poscounts path via +#' de_choose_size_factor_type) +#' - Both condition groups have identical counts (no differential signal, +#' stresses dispersion estimation) +#' +#' @param out_dir Directory to write the output files. +#' @return A list with paths to transcript RDS, gene RDS, and sample sheet CSV. +#' @export +write_degenerate_de_fixture_bundle <- function(out_dir) { + sample_names <- c( + "ctrl_rep1", "ctrl_rep2", "ctrl_rep3", + "case_rep1", "case_rep2", "case_rep3" + ) + + # Identical counts across both groups; zeros in every transcript. + tx_counts <- matrix( + c( + 0, 0, 4, 0, 0, 4, # tx1: zero in rep1 and rep2 of both groups + 3, 4, 3, 3, 4, 3, # tx2: all nonzero but no group difference + 2, 2, 0, 2, 2, 0, # tx3: zero in rep3 + 0, 1, 1, 0, 1, 1 # tx4: zero in rep1 + ), + nrow = 4, byrow = TRUE, + dimnames = list(c("tx1", "tx2", "tx3", "tx4"), sample_names) + ) + tx_cpm <- t(t(tx_counts) / pmax(colSums(tx_counts), 1)) * 1e6 + + row_ranges <- GenomicRanges::GRanges( + seqnames = rep("chr1", 4), + ranges = IRanges::IRanges(start = c(1, 101, 201, 301), width = 50), + strand = rep("+", 4), + TXNAME = c("tx1", "tx2", "tx3", "tx4"), + GENEID = c("gene1", "gene1", "gene2", "gene2"), + eqClassById = IRanges::CharacterList(list(c("1", "2"), "3", "4", "5")) + ) + + tx_se <- SummarizedExperiment::SummarizedExperiment( + assays = list(counts = tx_counts, CPM = tx_cpm), + rowRanges = row_ranges + ) + + # Gene counts: gene1 = tx1+tx2, gene2 = tx3+tx4. + # Each gene has at least one zero per group (triggers poscounts for DGE). + gene_counts <- rbind( + gene1 = tx_counts["tx1", ] + tx_counts["tx2", ], + gene2 = tx_counts["tx3", ] + tx_counts["tx4", ] + ) + gene_cpm <- t(t(gene_counts) / pmax(colSums(gene_counts), 1)) * 1e6 + + gene_se <- SummarizedExperiment::SummarizedExperiment( + assays = list(counts = gene_counts, CPM = gene_cpm), + rowData = S4Vectors::DataFrame(GENEID = rownames(gene_counts)) + ) + + tx_path <- file.path(out_dir, "transcripts.rds") + gene_path <- file.path(out_dir, "genes.rds") + saveRDS(tx_se, tx_path) + saveRDS(gene_se, gene_path) + + sample_sheet <- file.path(out_dir, "sample_sheet.csv") + utils::write.csv( + data.frame( + alias = sample_names, + condition = rep(c("control", "case"), each = 3), + stringsAsFactors = FALSE + ), + sample_sheet, + row.names = FALSE, + quote = FALSE + ) + + list( + transcript_rds = tx_path, + gene_rds = gene_path, + sample_sheet = sample_sheet + ) +} 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 ded16d2..b404eec 100644 --- a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R +++ b/bin/workflow_glue_r/tests/testthat/test_de_analysis.R @@ -535,9 +535,88 @@ testthat::test_that("de_is_recoverable_dexseq_error recognises expected messages )) 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")) + # Internal object-shape mismatches should still be reported as failures, + # but they are not known statistical fallback patterns. + testthat::expect_false(de_is_recoverable_dexseq_error( + "nb of cols in 'assay' (7) must equal nb of rows in 'colData' (12)" + )) + testthat::expect_true(de_is_recoverable_dexseq_error( + "newsplit: out of vertex space" + )) + testthat::expect_true(de_is_recoverable_dexseq_error( + "every gene contains at least one zero, cannot compute log geometric means" + )) testthat::expect_false(de_is_recoverable_dexseq_error("random unrelated failure")) }) +testthat::test_that("de_dexseq_failure_hint identifies degenerate-data shape mismatch", { + hint <- de_dexseq_failure_hint( + "nb of cols in 'assay' (7) must equal nb of rows in 'colData' (12)" + ) + testthat::expect_match(hint, "degenerate or near-identical sample groups") + row_sums_hint <- de_dexseq_failure_hint( + "Error in rowSums(numerator) : 'x' must be an array of at least two dimensions" + ) + testthat::expect_match(row_sums_hint, "gene-level DTU summarisation") + testthat::expect_null(de_dexseq_failure_hint("random unrelated failure")) +}) + +testthat::test_that("de_postprocess_dexseq_result degrades gracefully on gene-level DTU failure", { + contrast_dir <- tempfile("dexseq-postprocess-") + dir.create(contrast_dir) + + dex_res <- list( + dxr = data.frame( + featureID = c("tx1", "tx2"), + groupID = c("gene1", "gene1"), + log2fold_case_control = c(1, -1), + pvalue = c(0.01, 0.2), + padj = c(0.02, 0.3), + exonBaseMean = c(10, 20), + check.names = FALSE + ), + dxd = structure(list(), class = "mock_dxd") + ) + contrast_samples <- data.frame( + alias = c("s1", "s2"), + condition = c("case", "control"), + stringsAsFactors = FALSE + ) + tx_counts <- matrix( + c(10, 0, 4, 8), + nrow = 2, + dimnames = list(c("tx1", "tx2"), c("s1", "s2")) + ) + + result <- testthat::expect_warning( + de_postprocess_dexseq_result( + dex_res = dex_res, + contrast_dir = contrast_dir, + target_level = "case", + reference_level = "control", + contrast_samples = contrast_samples, + condition_column = "condition", + tx_counts = tx_counts, + contrast_name = "case_control", + per_gene_qvalue_fn = function(...) { + stop("'x' must be an array of at least two dimensions", call. = FALSE) + }, + plot_writer = function(...) NULL + ), + "DEXSeq post-processing failed" + ) + + testthat::expect_false(result$ok) + testthat::expect_match( + result$failure_hint, + "gene-level DTU summarisation" + ) + testthat::expect_true(file.exists(file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt"))) + failure_text <- readLines(file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt")) + testthat::expect_true(any(grepl("DTU Analysis Failed", failure_text, fixed = TRUE))) + testthat::expect_true(any(grepl("at least two dimensions", failure_text, fixed = TRUE))) +}) + testthat::test_that("de_run_dexseq_result records rank-deficiency covariate drops", { testthat::skip_if_not_installed("DESeq2") testthat::skip_if_not_installed("DEXSeq") @@ -779,7 +858,89 @@ testthat::test_that("CLI integration produces expected outputs", { contrast_qc <- de_qc$contrasts[["condition_treated_vs_control"]] testthat::expect_true("deseq2_size_factor_method" %in% names(contrast_qc)) testthat::expect_true("deseq2_dispersion_fallback" %in% names(contrast_qc)) + testthat::expect_true("dge_status" %in% names(contrast_qc)) testthat::expect_true("dexseq_size_factor_method" %in% names(contrast_qc)) testthat::expect_true("dexseq_dispersion_method" %in% names(contrast_qc)) testthat::expect_true("dexseq_covariates_dropped" %in% names(contrast_qc)) }) + +### +# Degenerate input (control-vs-control / zero-heavy counts) +# +# Exercises the graceful-degradation paths added for three cases: +# - newsplit/lfproc OOM in DESeq2 local regression +# - DEXSeq colData dimension mismatch on degenerate input +# - size-factor estimation failure when every gene has a zero +# +# The process must run to completion and produce placeholder outputs + a +# machine-readable status in de_qc_stats.json rather than crashing. + +testthat::test_that("degenerate control-vs-control data completes gracefully", { + testthat::skip_if_not_installed("DESeq2") + testthat::skip_if_not_installed("DEXSeq") + + fixture_dir <- tempfile("de-degenerate-") + dir.create(fixture_dir) + bundle <- write_degenerate_de_fixture_bundle(fixture_dir) + + argv <- c( + bundle, + list( + condition_column = "condition", + covariates = NULL, + reference_level = "control", + out_dir = file.path(fixture_dir, "out") + ) + ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) + + # Must complete without throwing — any DESeq2/DEXSeq failure should be + # caught and degraded to placeholder outputs rather than crashing the process. + result <- tryCatch( + suppressWarnings(suppressMessages(main_run_de_analysis(argv))), + error = function(err) err + ) + testthat::expect_false( + inherits(result, "error"), + info = if (inherits(result, "error")) conditionMessage(result) else "" + ) + + # QC JSON must exist with per-contrast status fields. + qc_path <- file.path(argv$out_dir, "de_qc_stats.json") + testthat::expect_true(file.exists(qc_path)) + de_qc <- jsonlite::read_json(qc_path, simplifyVector = TRUE) + + contrast_name <- "condition_case_vs_control" + contrast_qc <- de_qc$contrasts[[contrast_name]] + testthat::expect_false(is.null(contrast_qc)) + testthat::expect_true(contrast_qc$dge_status %in% c("SUCCESS", "FAILED")) + testthat::expect_true(contrast_qc$dtu_status %in% c("SUCCESS", "FAILED")) + + # DGE results file must always exist (empty placeholder on failure). + contrast_dir <- file.path(argv$out_dir, contrast_name) + dge_tsv <- file.path(contrast_dir, "results_dge.tsv") + testthat::expect_true(file.exists(dge_tsv)) + + if (identical(contrast_qc$dge_status, "FAILED")) { + # Failure log and placeholder PDF written. + testthat::expect_true(file.exists(file.path(contrast_dir, "DGE_ANALYSIS_FAILED.txt"))) + testthat::expect_true(file.exists(file.path(contrast_dir, "results_dge.pdf"))) + # Empty TSV must have the expected column headers. + dge_df <- utils::read.delim(dge_tsv, check.names = FALSE) + testthat::expect_equal(nrow(dge_df), 0) + testthat::expect_true(all(c("GENEID", "log2FoldChange", "padj") %in% names(dge_df))) + } + + # DTU outputs must always exist. + testthat::expect_true(file.exists(file.path(contrast_dir, "results_dtu_transcript.tsv"))) + testthat::expect_true(file.exists(file.path(contrast_dir, "results_dtu_gene.tsv"))) + + if (identical(contrast_qc$dtu_status, "FAILED")) { + testthat::expect_true(file.exists(file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt"))) + } + + # analysis_fallbacks must record failed counts. + fallbacks <- de_qc$analysis_fallbacks + testthat::expect_true(!is.null(fallbacks$failed_dge_contrasts)) + testthat::expect_true(!is.null(fallbacks$failed_dtu_contrasts)) +})