From 5af948dabf74142d2f26d341d63a06bca71d3987 Mon Sep 17 00:00:00 2001 From: Chris Wright Date: Thu, 14 May 2026 17:41:21 +0000 Subject: [PATCH] [CW-7214] Tidy up CLI parsing --- .gitlab-ci.yml | 14 +- bin/workflow_glue_r/NAMESPACE | 9 +- bin/workflow_glue_r/R/args.R | 213 +++++++++ bin/workflow_glue_r/R/bambu.R | 433 ++++++++---------- bin/workflow_glue_r/R/common.R | 39 -- bin/workflow_glue_r/R/de_analysis.R | 188 ++++---- .../tests/testthat/helper-fixtures.R | 65 +++ .../tests/testthat/test_bambu.R | 142 +++--- .../tests/testthat/test_de_analysis.R | 12 +- subworkflows/transcriptome.nf | 2 +- 10 files changed, 674 insertions(+), 443 deletions(-) create mode 100644 bin/workflow_glue_r/R/args.R diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index f88ffb8..c300708 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -57,7 +57,7 @@ docker-run: parallel: matrix: - MATRIX_NAME: [ - "de-poscounts-fallback", "discover", "igv", + "discover", "igv", "smoke_discover", "smoke_fixed", "smoke_direct_rna", "smoke_de", "no_annotation", "invalid_mode", "conflicting_flags" ] @@ -81,18 +81,6 @@ docker-run: -c ${CI_PROJECT_NAME}/data/demo.nextflow.config " ASSERT_NEXTFLOW_FAILURE: "1" ASSERT_NEXTFLOW_FAILURE_REXP: "Missing required parameter: --ref_annotation" - - if: $MATRIX_NAME == "de-poscounts-fallback" - variables: - NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && 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;" - 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 --minimap2_index_opts '-k 15' --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression/sample_sheet.csv \ - -c ${CI_PROJECT_NAME}/data/demo.nextflow.config " - NF_IGNORE_PROCESSES: faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome - AFTER_NEXTFLOW_CMD: > - grep -Eq '"deseq2_size_factor_method": "poscounts"' ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json - if: $MATRIX_NAME == "only_differential_expression" variables: NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && 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;" diff --git a/bin/workflow_glue_r/NAMESPACE b/bin/workflow_glue_r/NAMESPACE index 6d935f5..e22afa7 100644 --- a/bin/workflow_glue_r/NAMESPACE +++ b/bin/workflow_glue_r/NAMESPACE @@ -4,23 +4,18 @@ export(bambu_discovery_enabled) export(bambu_filter_transcripts) export(bambu_normalise_tsv_df) export(bambu_resolve_inputs) -export(bambu_resolve_ndr) export(bambu_strip_alias) -export(bambu_validate_args) export(bambu_write_outputs) export(de_analysis_arg_parser) -export(de_build_contrast_name) -export(de_parse_covariates) export(de_validate_inputs) export(main_run_bambu) export(main_run_de_analysis) export(run_bambu_cli) export(run_de_analysis_cli) +export(workflow_glue_r_arg_parser_from_spec) export(workflow_glue_r_cli) export(workflow_glue_r_components) -export(workflow_glue_r_arg_missing) export(workflow_glue_r_empty_tsv) +export(workflow_glue_r_normalise_args) export(workflow_glue_r_normalise_tsv_df) export(workflow_glue_r_parse_csv_list) -export(workflow_glue_r_read_csv) -export(workflow_glue_r_require_args) diff --git a/bin/workflow_glue_r/R/args.R b/bin/workflow_glue_r/R/args.R new file mode 100644 index 0000000..30b2b5c --- /dev/null +++ b/bin/workflow_glue_r/R/args.R @@ -0,0 +1,213 @@ +workflow_glue_r_parse_csv_list <- function(value) { + if (is.null(value) || length(value) == 0 || (length(value) == 1 && is.na(value))) { + return(character(0)) + } + + values <- trimws(strsplit(value, ",", fixed = TRUE)[[1]]) + values[nzchar(values)] +} + +workflow_glue_r_flag_present <- function(raw_argv, flag) { + if (is.null(raw_argv) || length(raw_argv) == 0) { + return(FALSE) + } + any(raw_argv == flag | startsWith(raw_argv, paste0(flag, "="))) +} + +workflow_glue_r_arg_parser_from_spec <- function(description, arg_spec) { + parser <- argparser::arg_parser(description) + for (arg in arg_spec) { + add_args <- list( + parser = parser, + arg = arg$flag, + help = arg$help, + type = arg$type + ) + if ("default" %in% names(arg)) { + add_args$default <- arg$default + } + parser <- do.call(argparser::add_argument, add_args) + } + parser +} + +workflow_glue_r_arg_value_error <- function(arg, fallback) { + if ("value_error" %in% names(arg)) { + return(arg$value_error) + } + fallback +} + +workflow_glue_r_scalar_arg <- function(value, arg) { + if (length(value) != 1) { + stop( + sprintf( + "%s must be a single %s value.", + arg$flag, + arg$type + ), + call. = FALSE + ) + } + value +} + +workflow_glue_r_normalise_arg_value <- function(value, arg, flag_provided = FALSE) { + value_is_na <- length(value) == 1 && is.na(value) + value_is_absent <- is.null(value) || + length(value) == 0 || + (value_is_na && !isTRUE(flag_provided)) || + (is.character(value) && length(value) == 1 && !nzchar(value)) + + if (value_is_absent) { + if ("default" %in% names(arg)) { + return(arg$default) + } + return(NULL) + } + + value <- workflow_glue_r_scalar_arg(value, arg) + if (identical(arg$type, "character")) { + if (is.na(value)) { + stop( + workflow_glue_r_arg_value_error(arg, sprintf("%s must be a non-empty string.", arg$flag)), + call. = FALSE + ) + } + return(as.character(value)) + } + + if (identical(arg$type, "integer")) { + numeric_value <- suppressWarnings(as.numeric(value)) + integer_value <- suppressWarnings(as.integer(numeric_value)) + if ( + is.na(numeric_value) || + is.na(integer_value) || + !is.finite(numeric_value) || + numeric_value != integer_value + ) { + stop( + workflow_glue_r_arg_value_error(arg, sprintf("%s must be an integer.", arg$flag)), + call. = FALSE + ) + } + return(integer_value) + } + + if (identical(arg$type, "numeric")) { + numeric_value <- suppressWarnings(as.numeric(value)) + if (is.na(numeric_value)) { + stop( + workflow_glue_r_arg_value_error(arg, sprintf("%s must be numeric.", arg$flag)), + call. = FALSE + ) + } + return(numeric_value) + } + + if (identical(arg$type, "logical")) { + logical_value <- suppressWarnings(as.logical(value)) + if (length(logical_value) != 1 || is.na(logical_value)) { + stop( + workflow_glue_r_arg_value_error(arg, sprintf("%s must be true or false.", arg$flag)), + call. = FALSE + ) + } + return(logical_value) + } + + value +} + +workflow_glue_r_validate_arg_value <- function(value, arg) { + if (is.null(value)) { + return(invisible(value)) + } + + if ("choices" %in% names(arg) && !value %in% arg$choices) { + stop( + sprintf( + "%s must be one of: %s", + arg$name, + paste(arg$choices, collapse = ", ") + ), + call. = FALSE + ) + } + + if ( + ("min" %in% names(arg) && value < arg$min) || + ("max" %in% names(arg) && value > arg$max) + ) { + stop( + workflow_glue_r_arg_value_error( + arg, + sprintf( + "%s must be between %s and %s.", + arg$flag, + if ("min" %in% names(arg)) arg$min else "-Inf", + if ("max" %in% names(arg)) arg$max else "Inf" + ) + ), + call. = FALSE + ) + } + + invisible(value) +} + +workflow_glue_r_normalise_args <- function(argv, arg_spec, raw_argv = NULL) { + for (arg in arg_spec) { + argv[[arg$name]] <- workflow_glue_r_normalise_arg_value( + argv[[arg$name]], + arg, + flag_provided = workflow_glue_r_flag_present(raw_argv, arg$flag) + ) + workflow_glue_r_validate_arg_value(argv[[arg$name]], arg) + } + + required_args <- vapply(arg_spec, function(arg) { + isTRUE(arg$required) + }, logical(1)) + required_arg_names <- vapply(arg_spec[required_args], `[[`, character(1), "name") + missing_args <- required_arg_names[vapply(required_arg_names, function(arg_name) { + is.null(argv[[arg_name]]) + }, logical(1))] + + if (length(missing_args) > 0) { + stop( + sprintf( + "Missing required arguments: %s", + paste(sprintf("--%s", missing_args), collapse = ", ") + ), + call. = FALSE + ) + } + + xor_args <- vapply(arg_spec, function(arg) { + "xor_group" %in% names(arg) + }, logical(1)) + xor_groups <- unique(vapply(arg_spec[xor_args], `[[`, character(1), "xor_group")) + for (xor_group in xor_groups) { + group_args <- arg_spec[vapply(arg_spec, function(arg) { + identical(arg$xor_group, xor_group) + }, logical(1))] + group_arg_names <- vapply(group_args, `[[`, character(1), "name") + present_args <- vapply(group_arg_names, function(arg_name) { + !is.null(argv[[arg_name]]) + }, logical(1)) + + if (sum(present_args) != 1) { + group_flags <- vapply(group_args, `[[`, character(1), "flag") + stop( + sprintf( + "Provide exactly one of %s.", + paste(group_flags, collapse = " or ") + ), + call. = FALSE + ) + } + } + + argv +} diff --git a/bin/workflow_glue_r/R/bambu.R b/bin/workflow_glue_r/R/bambu.R index 3555dc9..1e4c5b8 100644 --- a/bin/workflow_glue_r/R/bambu.R +++ b/bin/workflow_glue_r/R/bambu.R @@ -1,66 +1,89 @@ -bambu_arg_parser <- function() { - parser <- argparser::arg_parser("Run bambu transcript discovery and quantification.") - parser <- argparser::add_argument(parser, "--bams", help = "Comma-separated BAM paths.") - parser <- argparser::add_argument(parser, "--aliases", help = "Comma-separated aliases for --bams.") - parser <- argparser::add_argument(parser, "--sample_sheet", help = "Optional sample sheet CSV.") - parser <- argparser::add_argument(parser, "--annotation", help = "Reference annotation GTF/GFF.") - parser <- argparser::add_argument(parser, "--genome", help = "Reference genome FASTA.") - parser <- argparser::add_argument( - parser, - "--transcriptome_mode", - help = "discover or fixed_annotation.", - default = "discover" - ) - parser <- argparser::add_argument( - parser, - "--threads", - help = "Number of worker threads.", - type = "numeric", - default = 1 - ) - parser <- argparser::add_argument( - parser, - "--ndr", - help = "Optional novel discovery rate.", - type = "numeric" - ) - argparser::add_argument(parser, "--out_dir", help = "Output directory.") -} - -bambu_validate_args <- function(argv) { - workflow_glue_r_require_args(argv, c("annotation", "genome", "out_dir")) - - if (workflow_glue_r_arg_missing(argv$bams)) { - stop("Missing required arguments: --bams", call. = FALSE) - } - if (workflow_glue_r_arg_missing(argv$aliases)) { - stop("Missing required arguments: --aliases", call. = FALSE) - } - - if (!argv$transcriptome_mode %in% c("discover", "fixed_annotation")) { - stop( - sprintf( - "transcriptome_mode must be one of: %s", - paste(c("discover", "fixed_annotation"), collapse = ", ") - ), - call. = FALSE +bambu_arg_spec <- function() { + list( + list( + name = "bams", + flag = "--bams", + help = "Comma-separated BAM paths.", + type = "character", + required = TRUE + ), + list( + name = "aliases", + flag = "--aliases", + help = "Comma-separated aliases for --bams.", + type = "character", + required = TRUE + ), + list( + name = "sample_sheet", + flag = "--sample_sheet", + help = "Optional sample sheet CSV.", + type = "character" + ), + list( + name = "annotation", + flag = "--annotation", + help = "Reference annotation GTF/GFF.", + type = "character", + required = TRUE + ), + list( + name = "genome", + flag = "--genome", + help = "Reference genome FASTA.", + type = "character", + required = TRUE + ), + list( + name = "transcriptome_mode", + flag = "--transcriptome_mode", + help = "discover or fixed_annotation.", + type = "character", + default = "discover", + choices = c("discover", "fixed_annotation") + ), + list( + name = "threads", + flag = "--threads", + help = "Number of worker threads.", + type = "integer", + default = 1L, + min = 1L + ), + list( + name = "ndr", + flag = "--ndr", + help = "Optional novel discovery rate.", + type = "numeric", + min = 0, + max = 1, + value_error = "NDR (Novel Discovery Rate) must be between 0 and 1" + ), + list( + name = "out_dir", + flag = "--out_dir", + help = "Output directory.", + type = "character", + required = TRUE ) - } - - if (!workflow_glue_r_arg_missing(argv$ndr) && (argv$ndr < 0 || argv$ndr > 1)) { - stop("NDR (Novel Discovery Rate) must be between 0 and 1", call. = FALSE) - } - - invisible(argv) + ) } -bambu_resolve_inputs <- function( - argv, - bamfile_list_ctor = Rsamtools::BamFileList -) { +bambu_arg_parser <- function() { + workflow_glue_r_arg_parser_from_spec( + "Run bambu transcript discovery and quantification.", + bambu_arg_spec() + ) +} + +bambu_resolve_inputs <- function(args) { sample_df <- NULL - if (!workflow_glue_r_arg_missing(argv$sample_sheet)) { - sample_df <- workflow_glue_r_read_csv(argv$sample_sheet) + if (!is.null(args$sample_sheet)) { + sample_df <- utils::read.csv( + args$sample_sheet, + check.names = FALSE, + stringsAsFactors = FALSE + ) if (!"alias" %in% names(sample_df)) { stop("Sample sheet must contain an 'alias' column.", call. = FALSE) } @@ -76,8 +99,8 @@ bambu_resolve_inputs <- function( } } - bam_paths <- workflow_glue_r_parse_csv_list(argv$bams) - aliases <- workflow_glue_r_parse_csv_list(argv$aliases) + bam_paths <- workflow_glue_r_parse_csv_list(args$bams) + aliases <- workflow_glue_r_parse_csv_list(args$aliases) if (length(bam_paths) < 1) { stop("No BAM files were provided in --bams.", call. = FALSE) @@ -116,7 +139,7 @@ bambu_resolve_inputs <- function( reads <- if (length(bam_paths) == 1) { bam_paths } else { - bamfile_list_ctor(bam_paths, yieldSize = 1000000) + Rsamtools::BamFileList(bam_paths, yieldSize = 250000L) } list( @@ -127,34 +150,47 @@ bambu_resolve_inputs <- function( ) } -bambu_discovery_enabled <- function(argv) { - identical(argv$transcriptome_mode, "discover") +bambu_discovery_enabled <- function(args) { + identical(args$transcriptome_mode, "discover") } -bambu_resolve_ndr <- function(argv, default_ndr = 0.1) { - if (workflow_glue_r_arg_missing(argv$ndr)) { - default_ndr - } else { - as.numeric(argv$ndr) - } -} - -bambu_build_args <- function(argv, reads, annotation_obj) { +bambu_build_args <- function(args, reads, annotation_obj) { bambu_args <- list( reads = reads, annotations = annotation_obj, - genome = argv$genome, - ncore = as.integer(argv$threads), - discovery = bambu_discovery_enabled(argv) + genome = args$genome, + ncore = args$threads, + discovery = bambu_discovery_enabled(args), + lowMemory = TRUE, + yieldSize = 250000L, + verbose = TRUE ) - if (bambu_discovery_enabled(argv)) { - bambu_args$NDR <- bambu_resolve_ndr(argv) + if (bambu_discovery_enabled(args) && !is.null(args$ndr)) { + bambu_args$NDR <- args$ndr } bambu_args } +bambu_effective_threads <- function(args, bam_count) { + # bambu's low-memory mode can have issues with multiple BAMs and + # parallel threads due to BiocFileCache writes, + # so we enforce single-threading in that case. + threads <- as.integer(args$threads) + if (bam_count > 1 && threads > 1L) { + warning( + paste( + "Low-memory mode with multiple BAMs can fail in bambu due to", + "parallel BiocFileCache writes; forcing threads=1." + ), + call. = FALSE + ) + return(1L) + } + threads +} + bambu_filter_transcripts <- function(se) { counts_mat <- SummarizedExperiment::assays(se)$counts full_length_mat <- SummarizedExperiment::assays(se)$fullLengthCounts @@ -194,85 +230,6 @@ bambu_matrix_to_df <- function(se_obj, assay_name, id_col, meta_df) { merge(meta_df, assay_df, by.x = id_col, by.y = id_col, all.y = TRUE, sort = FALSE) } -bambu_extract_gtf_attribute <- function(attr_field, key) { - match <- regexec(sprintf('%s "([^"]*)";', key), attr_field, perl = TRUE) - captures <- regmatches(attr_field, match)[[1]] - if (length(captures) < 2) { - return(NULL) - } - captures[2] -} - -bambu_normalise_gtf_attribute_value <- function(value) { - if (is.null(value)) { - return(NULL) - } - value <- gsub('[";]', "", value) - value <- trimws(gsub("\\s+", " ", value)) - if (!nzchar(value)) { - return(NULL) - } - value -} - -bambu_sanitise_gtf_file <- function(path) { - lines <- readLines(path, warn = FALSE) - cleaned_lines <- vapply(lines, function(line) { - if (!nzchar(line) || startsWith(line, "#")) { - return(line) - } - - fields <- strsplit(line, "\t", fixed = TRUE)[[1]] - if (length(fields) < 9) { - return(line) - } - - attr_field <- fields[9] - transcript_id <- bambu_normalise_gtf_attribute_value( - bambu_extract_gtf_attribute(attr_field, "transcript_id") - ) - gene_id <- bambu_normalise_gtf_attribute_value( - bambu_extract_gtf_attribute(attr_field, "gene_id") - ) - - if (!is.null(gene_id) && identical(gene_id, "transcript_id")) { - warning( - sprintf( - "Replaced malformed gene_id 'transcript_id' with transcript_id '%s'.", - transcript_id - ), - call. = FALSE - ) - gene_id <- transcript_id - } - if (is.null(gene_id)) { - gene_id <- transcript_id - } - - if (!is.null(gene_id)) { - attr_field <- sub( - 'gene_id "([^"]*)";', - sprintf('gene_id "%s";', gene_id), - attr_field, - perl = TRUE - ) - } - if (!is.null(transcript_id)) { - attr_field <- sub( - 'transcript_id "([^"]*)";', - sprintf('transcript_id "%s";', transcript_id), - attr_field, - perl = TRUE - ) - } - - fields[9] <- attr_field - paste(fields, collapse = "\t") - }, character(1)) - - writeLines(cleaned_lines, path) -} - bambu_format_count <- function(value) { if (length(value) == 0 || all(is.na(value))) { return("NA") @@ -285,17 +242,30 @@ bambu_format_count <- function(value) { ) } -bambu_write_outputs <- function(se, gene_se, sample_df, argv, qc_stats, write_gtf_fn = bambu::writeToGTF) { - write_gtf_fn( +bambu_write_matrix_tsv <- function(se_obj, assay_name, id_col, meta_df, output_path) { + table_df <- bambu_matrix_to_df(se_obj, assay_name, id_col, meta_df) + utils::write.table( + table_df, + file = output_path, + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + rm(table_df) + invisible(gc(verbose = FALSE)) +} + +bambu_write_outputs <- function(se, gene_se, sample_df, args, qc_stats) { + bambu::writeToGTF( SummarizedExperiment::rowRanges(se), - file = file.path(argv$out_dir, "transcripts.gtf") + file = file.path(args$out_dir, "transcripts.gtf") ) - saveRDS(se, file.path(argv$out_dir, "bambu_transcripts.rds")) - saveRDS(gene_se, file.path(argv$out_dir, "bambu_genes.rds")) + saveRDS(se, file.path(args$out_dir, "bambu_transcripts.rds")) + saveRDS(gene_se, file.path(args$out_dir, "bambu_genes.rds")) utils::write.csv( sample_df, - file.path(argv$out_dir, "samples.csv"), + file.path(args$out_dir, "samples.csv"), row.names = FALSE, quote = FALSE ) @@ -318,58 +288,51 @@ bambu_write_outputs <- function(se, gene_se, sample_df, argv, qc_stats, write_gt utils::write.table( tx_meta, - file = file.path(argv$out_dir, "transcript_metadata.tsv"), + file = file.path(args$out_dir, "transcript_metadata.tsv"), sep = "\t", quote = FALSE, row.names = FALSE ) utils::write.table( gene_meta, - file = file.path(argv$out_dir, "gene_metadata.tsv"), + file = file.path(args$out_dir, "gene_metadata.tsv"), sep = "\t", quote = FALSE, row.names = FALSE ) - tx_counts <- bambu_matrix_to_df(se, "counts", "TXNAME", tx_meta) - tx_cpm <- bambu_matrix_to_df(se, "CPM", "TXNAME", tx_meta) - gene_counts <- bambu_matrix_to_df(gene_se, "counts", "GENEID", gene_meta) - gene_cpm <- bambu_matrix_to_df(gene_se, "CPM", "GENEID", gene_meta) - - utils::write.table( - tx_counts, - file = file.path(argv$out_dir, "transcript_counts.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE + bambu_write_matrix_tsv( + se, + "counts", + "TXNAME", + tx_meta, + file.path(args$out_dir, "transcript_counts.tsv") ) - utils::write.table( - tx_cpm, - file = file.path(argv$out_dir, "transcript_cpm.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE + bambu_write_matrix_tsv( + se, + "CPM", + "TXNAME", + tx_meta, + file.path(args$out_dir, "transcript_cpm.tsv") ) - utils::write.table( - gene_counts, - file = file.path(argv$out_dir, "gene_counts.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE + bambu_write_matrix_tsv( + gene_se, + "counts", + "GENEID", + gene_meta, + file.path(args$out_dir, "gene_counts.tsv") ) - utils::write.table( - gene_cpm, - file = file.path(argv$out_dir, "gene_cpm.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE + bambu_write_matrix_tsv( + gene_se, + "CPM", + "GENEID", + gene_meta, + file.path(args$out_dir, "gene_cpm.tsv") ) - bambu_sanitise_gtf_file(file.path(argv$out_dir, "transcripts.gtf")) - - qc_stats$transcriptome_mode <- argv$transcriptome_mode - qc_stats$ndr_used <- if (bambu_discovery_enabled(argv)) { - bambu_resolve_ndr(argv) + qc_stats$transcriptome_mode <- args$transcriptome_mode + qc_stats$ndr_used <- if (bambu_discovery_enabled(args)) { + if (is.null(args$ndr)) "automatic" else args$ndr } else { "N/A" } @@ -377,7 +340,7 @@ bambu_write_outputs <- function(se, gene_se, sample_df, argv, qc_stats, write_gt jsonlite::write_json( qc_stats, - file.path(argv$out_dir, "bambu_qc_stats.json"), + file.path(args$out_dir, "bambu_qc_stats.json"), pretty = TRUE, auto_unbox = TRUE ) @@ -387,8 +350,14 @@ bambu_write_outputs <- function(se, gene_se, sample_df, argv, qc_stats, write_gt "================================", "", sprintf("Timestamp: %s", qc_stats$timestamp), - sprintf("Mode: %s", argv$transcriptome_mode), - if (bambu_discovery_enabled(argv)) sprintf("NDR: %.3f", bambu_resolve_ndr(argv)) else NULL, + sprintf("Mode: %s", args$transcriptome_mode), + if (bambu_discovery_enabled(args)) { + if (is.null(args$ndr)) { + "NDR: automatic (bambu-selected)" + } else { + sprintf("NDR: %.3f", args$ndr) + } + } else NULL, "", "Sample Statistics:", sprintf(" Samples analyzed: %s", bambu_format_count(qc_stats$samples)), @@ -415,49 +384,51 @@ bambu_write_outputs <- function(se, gene_se, sample_df, argv, qc_stats, write_gt "" ) - writeLines(qc_summary, file.path(argv$out_dir, "bambu_qc_summary.txt")) - writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt")) + writeLines(qc_summary, file.path(args$out_dir, "bambu_qc_summary.txt")) + writeLines(capture.output(sessionInfo()), file.path(args$out_dir, "session_info.txt")) } main_run_bambu <- function( - argv, + args, analysis_fn = bambu::bambu, prepare_annotations_fn = bambu::prepareAnnotations, - gene_expression_fn = bambu::transcriptToGeneExpression, - write_gtf_fn = bambu::writeToGTF, - bamfile_list_ctor = Rsamtools::BamFileList + gene_expression_fn = bambu::transcriptToGeneExpression ) { set.seed(42) - suppressPackageStartupMessages({ - library(GenomicRanges) - library(Rsamtools) - }) + # bambu's parallel worker code may rely on these being attached for generics + # such as seqlengths(). + suppressPackageStartupMessages(library(GenomicRanges)) + suppressPackageStartupMessages(library(Rsamtools)) - bambu_validate_args(argv) - dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE) + dir.create(args$out_dir, showWarnings = FALSE, recursive = TRUE) - inputs <- bambu_resolve_inputs( - argv, - bamfile_list_ctor = bamfile_list_ctor - ) - annotation_obj <- prepare_annotations_fn(argv$annotation) - ndr_value <- bambu_resolve_ndr(argv) + inputs <- bambu_resolve_inputs(args) + args$threads <- bambu_effective_threads(args, length(inputs$bam_paths)) + annotation_obj <- prepare_annotations_fn(args$annotation) - if (!workflow_glue_r_arg_missing(argv$ndr)) { - message(sprintf("Using user-specified NDR = %.3f", ndr_value)) - } else { - message(sprintf("Using default NDR = %.3f", ndr_value)) + if (!is.null(args$ndr)) { + message(sprintf("Using user-specified NDR = %.3f", args$ndr)) + } else if (bambu_discovery_enabled(args)) { + message("Using bambu automatic NDR selection.") } - if (bambu_discovery_enabled(argv)) { + if (bambu_discovery_enabled(args)) { message("Novel Discovery Rate (NDR) controls transcript discovery stringency:") message(" Lower NDR (e.g., 0.05) = fewer false positive transcripts, may miss real ones") message(" Higher NDR (e.g., 0.2) = more sensitive discovery, more false positives") - message(sprintf(" Current NDR = %.3f balances precision and recall", ndr_value)) + if (is.null(args$ndr)) { + message(" Current NDR = automatic (selected by bambu from the data)") + } else { + message(sprintf(" Current NDR = %.3f balances precision and recall", args$ndr)) + } } + if (length(inputs$bam_paths) > 1) { + message("Using BamFileList yieldSize = 250000") + } + message(sprintf("Running bambu with threads = %d", args$threads)) message("Running bambu...") - se <- do.call(analysis_fn, bambu_build_args(argv, inputs$reads, annotation_obj)) + se <- do.call(analysis_fn, bambu_build_args(args, inputs$reads, annotation_obj)) message("Bambu completed successfully") colnames(se) <- inputs$aliases @@ -509,9 +480,8 @@ main_run_bambu <- function( se, gene_se, inputs$sample_df, - argv, - qc_stats, - write_gtf_fn = write_gtf_fn + args, + qc_stats ) invisible( @@ -526,5 +496,6 @@ main_run_bambu <- function( run_bambu_cli <- function(argv = commandArgs(trailingOnly = TRUE)) { parsed <- argparser::parse_args(bambu_arg_parser(), argv = argv) - main_run_bambu(parsed) + args <- workflow_glue_r_normalise_args(parsed, bambu_arg_spec(), raw_argv = argv) + main_run_bambu(args) } diff --git a/bin/workflow_glue_r/R/common.R b/bin/workflow_glue_r/R/common.R index 0fa0caa..c78524e 100644 --- a/bin/workflow_glue_r/R/common.R +++ b/bin/workflow_glue_r/R/common.R @@ -1,38 +1,3 @@ -workflow_glue_r_arg_missing <- function(value) { - if (is.null(value) || length(value) == 0 || all(is.na(value))) { - return(TRUE) - } - if (is.character(value)) { - return(all(!nzchar(value))) - } - FALSE -} - -workflow_glue_r_require_args <- function(argv, required_args) { - missing_args <- required_args[vapply(required_args, function(arg_name) { - workflow_glue_r_arg_missing(argv[[arg_name]]) - }, logical(1))] - - if (length(missing_args) > 0) { - stop( - sprintf( - "Missing required arguments: %s", - paste(sprintf("--%s", missing_args), collapse = ", ") - ), - call. = FALSE - ) - } -} - -workflow_glue_r_parse_csv_list <- function(value) { - if (workflow_glue_r_arg_missing(value)) { - return(character(0)) - } - - values <- trimws(strsplit(value, ",", fixed = TRUE)[[1]]) - values[nzchar(values)] -} - workflow_glue_r_is_r_formula_name <- function(name) { is.character(name) && length(name) == 1 && @@ -54,10 +19,6 @@ workflow_glue_r_validate_r_formula_names <- function(names, label = "Column") { invisible(names) } -workflow_glue_r_read_csv <- function(path) { - utils::read.csv(path, check.names = FALSE, stringsAsFactors = FALSE) -} - workflow_glue_r_normalise_tsv_value <- function(value) { if (length(value) == 0 || all(is.na(value))) { return(NA_character_) diff --git a/bin/workflow_glue_r/R/de_analysis.R b/bin/workflow_glue_r/R/de_analysis.R index b096cdf..0205f92 100644 --- a/bin/workflow_glue_r/R/de_analysis.R +++ b/bin/workflow_glue_r/R/de_analysis.R @@ -1,27 +1,64 @@ -de_analysis_arg_parser <- function() { - parser <- argparser::arg_parser("Run DESeq2 and DEXSeq on bambu output.") - parser <- argparser::add_argument(parser, "--transcript_rds", help = "bambu transcript RDS.") - parser <- argparser::add_argument(parser, "--gene_rds", help = "bambu gene RDS.") - parser <- argparser::add_argument(parser, "--sample_sheet", help = "Sample sheet CSV.") - parser <- argparser::add_argument( - parser, - "--condition_column", - help = "Primary condition column.", - default = "condition" +de_analysis_arg_spec <- function() { + list( + list( + name = "transcript_rds", + flag = "--transcript_rds", + help = "bambu transcript RDS.", + type = "character", + required = TRUE + ), + list( + name = "gene_rds", + flag = "--gene_rds", + help = "bambu gene RDS.", + type = "character", + required = TRUE + ), + list( + name = "sample_sheet", + flag = "--sample_sheet", + help = "Sample sheet CSV.", + type = "character", + required = TRUE + ), + list( + name = "condition_column", + flag = "--condition_column", + help = "Primary condition column.", + type = "character", + default = "condition" + ), + list( + name = "covariates", + flag = "--covariates", + help = "Comma-separated nuisance covariates.", + type = "character" + ), + list( + name = "reference_level", + flag = "--reference_level", + help = "Reference level for the condition column.", + type = "character" + ), + list( + name = "out_dir", + flag = "--out_dir", + help = "Output directory.", + type = "character", + default = "de_analysis" + ) ) - parser <- argparser::add_argument(parser, "--covariates", help = "Comma-separated nuisance covariates.") - parser <- argparser::add_argument(parser, "--reference_level", help = "Reference level for the condition column.") - argparser::add_argument(parser, "--out_dir", help = "Output directory.", default = "de_analysis") } -de_parse_covariates <- function(value) { - workflow_glue_r_parse_csv_list(value) +de_analysis_arg_parser <- function() { + workflow_glue_r_arg_parser_from_spec( + "Run DESeq2 and DEXSeq on bambu output.", + de_analysis_arg_spec() + ) } de_validate_inputs <- function(tx_se, gene_se, sample_df, argv) { - workflow_glue_r_require_args(argv, c("transcript_rds", "gene_rds", "sample_sheet")) - - covariates <- de_parse_covariates(argv$covariates) + covariates <- workflow_glue_r_parse_csv_list(argv$covariates) workflow_glue_r_validate_r_formula_names( c(argv$condition_column, covariates), label = "Design column" @@ -107,7 +144,7 @@ de_validate_inputs <- function(tx_se, gene_se, sample_df, argv) { } reference_level <- argv$reference_level - if (workflow_glue_r_arg_missing(reference_level)) { + if (is.null(reference_level)) { if ("control" %in% condition_values) { reference_level <- "control" } else { @@ -127,8 +164,6 @@ de_validate_inputs <- function(tx_se, gene_se, sample_df, argv) { } list( - tx_se = tx_se, - gene_se = gene_se, sample_df = sample_df, covariates = covariates, condition_values = condition_values, @@ -136,19 +171,6 @@ de_validate_inputs <- function(tx_se, gene_se, sample_df, argv) { ) } -de_build_contrast_name <- function(condition_column, target_level, reference_level) { - sprintf("%s_%s_vs_%s", condition_column, target_level, reference_level) -} - -de_set_dispersions <- function(object, value) { - setter <- get("dispersions<-", envir = asNamespace("DESeq2")) - setter(object, value = value) -} - -de_extract_disp_gene_est <- function(object) { - S4Vectors::mcols(object)$dispGeneEst -} - # DESeq2's default geometric-mean size-factor estimator is undefined when # every gene has at least one zero across samples, so use poscounts then. de_choose_size_factor_type <- function(count_mat, context_label = "Count matrix") { @@ -161,7 +183,6 @@ de_choose_size_factor_type <- function(count_mat, context_label = "Count matrix" } "ratio" } - de_run_deseq_with_fallback <- function( dds, contrast_name, @@ -204,7 +225,8 @@ de_run_deseq_with_fallback <- function( dds <- DESeq2::estimateSizeFactors(dds, type = sf_type) dds <- DESeq2::estimateDispersionsGeneEst(dds) - dds <- de_set_dispersions(dds, de_extract_disp_gene_est(dds)) + dispersions_setter <- get("dispersions<-", envir = asNamespace("DESeq2")) + dds <- dispersions_setter(dds, value = S4Vectors::mcols(dds)$dispGeneEst) dispersion_values <- suppressWarnings(as.numeric(DESeq2::dispersions(dds))) dispersion_values <- dispersion_values[is.finite(dispersion_values)] @@ -279,8 +301,7 @@ de_estimate_dispersions_with_fallback <- function( list( object = DESeq2::estimateDispersions(object), method_used = "parametric", - fallback_applied = FALSE, - reason = NULL + fallback_applied = FALSE ), error = function(err) { if (!grepl( @@ -291,7 +312,6 @@ de_estimate_dispersions_with_fallback <- function( stop(err) } - primary_reason <- conditionMessage(err) message( context_label, " dispersion fitting failed; retrying with fitType='local'." @@ -300,8 +320,7 @@ de_estimate_dispersions_with_fallback <- function( list( object = DESeq2::estimateDispersions(object, fitType = "local"), method_used = "local", - fallback_applied = TRUE, - reason = primary_reason + fallback_applied = TRUE ), error = function(local_err) { if (!grepl( @@ -320,8 +339,7 @@ de_estimate_dispersions_with_fallback <- function( list( object = DESeq2::estimateDispersions(object, fitType = "mean"), method_used = "mean", - fallback_applied = TRUE, - reason = primary_reason + fallback_applied = TRUE ), error = function(mean_err) { if (!grepl( @@ -340,12 +358,15 @@ de_estimate_dispersions_with_fallback <- function( " mean-fit dispersion retry failed; falling back to gene-wise dispersion estimates." ) object <- DESeq2::estimateDispersionsGeneEst(object) - object <- de_set_dispersions(object, de_extract_disp_gene_est(object)) + dispersions_setter <- get("dispersions<-", envir = asNamespace("DESeq2")) + object <- dispersions_setter( + object, + value = S4Vectors::mcols(object)$dispGeneEst + ) list( object = object, method_used = "gene-wise", - fallback_applied = TRUE, - reason = primary_reason + fallback_applied = TRUE ) } ) @@ -404,7 +425,6 @@ de_run_deseq2_result <- function( independentFiltering = TRUE ) list( - dds = dds, result = result, deseq2_dispersion_fallback = deseq2_dispersion_fallback ) @@ -469,15 +489,8 @@ de_run_dexseq_result <- function( "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 <- dispersion_result$object + dispersion_method <- dispersion_result$method_used dxd <- DEXSeq::testForDEU(dxd, reducedModel = reduced_formula) dxd <- DEXSeq::estimateExonFoldChanges(dxd, fitExpToVar = condition_column) dxr <- DEXSeq::DEXSeqResults(dxd, independentFiltering = FALSE) @@ -485,7 +498,6 @@ de_run_dexseq_result <- function( dxd = dxd, dxr = dxr, dexseq_dispersion_method = dispersion_method, - dexseq_dispersion_reason = dispersion_reason, dexseq_size_factor_type = dexseq_sf_type ) }, error = function(err) { @@ -527,7 +539,7 @@ de_dtu_transcript_columns <- c( #' #' Renames the contrast-specific DEXSeq fold-change column to #' `log2FoldChange`, normalizes data for TSV output -#' and returns only the transcript output columns. +#' and returns only the transcript output columns. #' #' @param dex_df DEXSeq results as a data frame. #' @param contrast_name Contrast suffix used in the DEXSeq fold-change column. @@ -545,14 +557,18 @@ de_extract_dtu_transcript_table <- function(dex_df, contrast_name) { workflow_glue_r_normalise_tsv_df(tx_dtu) } -main_run_de_analysis <- function(argv) { +main_run_de_analysis <- function(args) { set.seed(42) - dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE) + dir.create(args$out_dir, showWarnings = FALSE, recursive = TRUE) - tx_se <- readRDS(argv$transcript_rds) - gene_se <- readRDS(argv$gene_rds) - sample_df <- workflow_glue_r_read_csv(argv$sample_sheet) - validated <- de_validate_inputs(tx_se, gene_se, sample_df, argv) + tx_se <- readRDS(args$transcript_rds) + gene_se <- readRDS(args$gene_rds) + sample_df <- utils::read.csv( + args$sample_sheet, + check.names = FALSE, + stringsAsFactors = FALSE + ) + validated <- de_validate_inputs(tx_se, gene_se, sample_df, args) sample_df <- validated$sample_df covariates <- validated$covariates condition_values <- validated$condition_values @@ -576,14 +592,14 @@ main_run_de_analysis <- function(argv) { de_qc_stats <- list( timestamp = format(Sys.time(), "%Y-%m-%d %H:%M:%S"), total_samples = nrow(sample_df), - condition_column = argv$condition_column, + condition_column = args$condition_column, reference_level = reference_level, covariates = if (length(covariates) > 0) covariates else "none", num_contrasts = length(targets), contrasts = list() ) - n_per_group <- table(sample_df[[argv$condition_column]]) + n_per_group <- table(sample_df[[args$condition_column]]) de_qc_stats$samples_per_group <- as.list(n_per_group) sample_size_warnings <- character(0) @@ -641,18 +657,23 @@ main_run_de_analysis <- function(argv) { " 4. Consider using hierarchical testing procedures", "" ) - writeLines(mt_content, file.path(argv$out_dir, "MULTIPLE_TESTING_WARNING.txt")) + writeLines(mt_content, file.path(args$out_dir, "MULTIPLE_TESTING_WARNING.txt")) } for (target_level in targets) { - contrast_name <- de_build_contrast_name(argv$condition_column, target_level, reference_level) - contrast_dir <- file.path(argv$out_dir, contrast_name) + contrast_name <- sprintf( + "%s_%s_vs_%s", + args$condition_column, + target_level, + reference_level + ) + contrast_dir <- file.path(args$out_dir, contrast_name) dir.create(contrast_dir, showWarnings = FALSE, recursive = TRUE) - keep_samples <- sample_df[[argv$condition_column]] %in% c(reference_level, target_level) + keep_samples <- sample_df[[args$condition_column]] %in% c(reference_level, target_level) contrast_samples <- droplevels(sample_df[keep_samples, , drop = FALSE]) - contrast_samples[[argv$condition_column]] <- stats::relevel( - factor(contrast_samples[[argv$condition_column]]), + contrast_samples[[args$condition_column]] <- stats::relevel( + factor(contrast_samples[[args$condition_column]]), ref = reference_level ) @@ -661,8 +682,8 @@ main_run_de_analysis <- function(argv) { target_level = target_level, 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_target = sum(contrast_samples[[args$condition_column]] == target_level), + n_reference = sum(contrast_samples[[args$condition_column]] == reference_level), deseq2_size_factor_method = "ratio", deseq2_dispersion_fallback = list( applied = FALSE, @@ -693,9 +714,9 @@ main_run_de_analysis <- function(argv) { contrast_samples, target_level, reference_level, - argv$condition_column, + args$condition_column, covariates, - argv$out_dir, + args$out_dir, contrast_name ) if (!is.null(dge_run$deseq2_dispersion_fallback)) { @@ -749,7 +770,7 @@ main_run_de_analysis <- function(argv) { tx_counts, tx_meta, contrast_samples, - argv$condition_column, + args$condition_column, covariates ), error = function(err) { @@ -776,9 +797,9 @@ main_run_de_analysis <- function(argv) { sprintf( "Samples: %d (%d %s, %d %s)", nrow(contrast_samples), - sum(contrast_samples[[argv$condition_column]] == target_level), + sum(contrast_samples[[args$condition_column]] == target_level), target_level, - sum(contrast_samples[[argv$condition_column]] == reference_level), + sum(contrast_samples[[args$condition_column]] == reference_level), reference_level ), sprintf("Transcripts: %d", nrow(tx_counts)), @@ -963,7 +984,7 @@ main_run_de_analysis <- function(argv) { jsonlite::write_json( de_qc_stats, - file.path(argv$out_dir, "de_qc_stats.json"), + file.path(args$out_dir, "de_qc_stats.json"), pretty = TRUE, auto_unbox = TRUE ) @@ -1008,13 +1029,14 @@ main_run_de_analysis <- function(argv) { " - /results_dtu_gene.tsv", "" ) - writeLines(unlist(overall_summary), file.path(argv$out_dir, "de_overall_summary.txt")) - writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt")) + writeLines(unlist(overall_summary), file.path(args$out_dir, "de_overall_summary.txt")) + writeLines(capture.output(sessionInfo()), file.path(args$out_dir, "session_info.txt")) invisible(list(qc = de_qc_stats)) } run_de_analysis_cli <- function(argv = commandArgs(trailingOnly = TRUE)) { parsed <- argparser::parse_args(de_analysis_arg_parser(), argv = argv) - main_run_de_analysis(parsed) + args <- workflow_glue_r_normalise_args(parsed, de_analysis_arg_spec(), raw_argv = argv) + main_run_de_analysis(args) } diff --git a/bin/workflow_glue_r/tests/testthat/helper-fixtures.R b/bin/workflow_glue_r/tests/testthat/helper-fixtures.R index 42307e1..eedfb08 100644 --- a/bin/workflow_glue_r/tests/testthat/helper-fixtures.R +++ b/bin/workflow_glue_r/tests/testthat/helper-fixtures.R @@ -118,6 +118,71 @@ make_test_tx_se <- function(include_geneid = TRUE, sample_names = NULL) { ) } +#' Create bambu-style transcript row ranges for output-writing tests. +#' +#' `bambu::writeToGTF()` expects a `GRangesList` shaped like the object +#' returned by `bambu::prepareAnnotations()`. This helper writes a minimal GTF +#' and returns that object with metadata columns used by test assertions. +#' +#' @param out_dir Directory where temporary annotation fixture is written. +#' @return A GRangesList with tx1-tx4 transcript entries and metadata. +#' @export +make_test_bambu_row_ranges <- function(out_dir) { + gtf <- file.path(out_dir, "annotation.gtf") + writeLines( + c( + paste( + "chr1", "test", "transcript", "1", "50", ".", "+", ".", + 'gene_id "gene1"; transcript_id "tx1";', + sep = "\t" + ), + paste( + "chr1", "test", "exon", "1", "50", ".", "+", ".", + 'gene_id "gene1"; transcript_id "tx1"; exon_number "1";', + sep = "\t" + ), + paste( + "chr1", "test", "transcript", "101", "150", ".", "+", ".", + 'gene_id "gene1"; transcript_id "tx2";', + sep = "\t" + ), + paste( + "chr1", "test", "exon", "101", "150", ".", "+", ".", + 'gene_id "gene1"; transcript_id "tx2"; exon_number "1";', + sep = "\t" + ), + paste( + "chr1", "test", "transcript", "201", "250", ".", "+", ".", + 'gene_id "gene2"; transcript_id "tx3";', + sep = "\t" + ), + paste( + "chr1", "test", "exon", "201", "250", ".", "+", ".", + 'gene_id "gene2"; transcript_id "tx3"; exon_number "1";', + sep = "\t" + ), + paste( + "chr1", "test", "transcript", "301", "350", ".", "+", ".", + 'gene_id "gene2"; transcript_id "tx4";', + sep = "\t" + ), + paste( + "chr1", "test", "exon", "301", "350", ".", "+", ".", + 'gene_id "gene2"; transcript_id "tx4"; exon_number "1";', + sep = "\t" + ) + ), + gtf + ) + + row_ranges <- bambu::prepareAnnotations(gtf)[c("tx1", "tx2", "tx3", "tx4")] + S4Vectors::mcols(row_ranges)$TXNAME <- names(row_ranges) + S4Vectors::mcols(row_ranges)$GENEID <- c("gene1", "gene1", "gene2", "gene2") + S4Vectors::mcols(row_ranges)$eqClassById <- IRanges::CharacterList(list(c("1", "2"), "3", "4", "5")) + + row_ranges +} + #' Create a test gene-level SummarizedExperiment by aggregating a transcript-level SummarizedExperiment. #' @param sample_names Optional vector of sample names to use as column names. #' @return A SummarizedExperiment object with synthetic gene-level data. diff --git a/bin/workflow_glue_r/tests/testthat/test_bambu.R b/bin/workflow_glue_r/tests/testthat/test_bambu.R index 2bd1140..17b2d6a 100644 --- a/bin/workflow_glue_r/tests/testthat/test_bambu.R +++ b/bin/workflow_glue_r/tests/testthat/test_bambu.R @@ -1,6 +1,6 @@ #' These tests cover the validation logic owned by supeRglue bambu before bambu -#' itself is invoked: mutually exclusive BAM inputs, alias derivation, sample -#' sheet alignment, transcriptome mode selection, and NDR handling. +#' itself is invoked: BAM/alias argument checks, sample-sheet alignment, +#' transcriptome mode selection, and NDR handling. #' #' NOTE: Annotation/reference preparation is handled by Python #' (bin/workflow_glue/prepare_annotation_reference.py) with pytest coverage. @@ -20,18 +20,20 @@ testthat::test_that("BAM inputs required", { ) testthat::expect_error( - bambu_validate_args(args), + workflow_glue_r_normalise_args(args, bambu_arg_spec()), "Missing required arguments: --bams" ) args$bams <- "sampleA.bam" testthat::expect_error( - bambu_validate_args(args), + workflow_glue_r_normalise_args(args, bambu_arg_spec()), "Missing required arguments: --aliases" ) args$aliases <- "sampleA" - testthat::expect_silent(bambu_validate_args(args)) + normalised <- NULL + testthat::expect_silent(normalised <- workflow_glue_r_normalise_args(args, bambu_arg_spec())) + testthat::expect_identical(normalised$threads, 1L) }) # transcriptome_mode must be "discover" or "fixed_annotation". @@ -48,27 +50,32 @@ testthat::test_that("invalid discovery settings rejected", { ) testthat::expect_error( - bambu_validate_args(args), + workflow_glue_r_normalise_args(args, bambu_arg_spec()), "transcriptome_mode must be one of" ) args$transcriptome_mode <- "discover" args$ndr <- -0.01 testthat::expect_error( - bambu_validate_args(args), + workflow_glue_r_normalise_args(args, bambu_arg_spec()), "NDR .* must be between 0 and 1" ) args$ndr <- 1.01 testthat::expect_error( - bambu_validate_args(args), + workflow_glue_r_normalise_args(args, bambu_arg_spec()), "NDR .* must be between 0 and 1" ) args$ndr <- 0 - testthat::expect_silent(bambu_validate_args(args)) + testthat::expect_silent(workflow_glue_r_normalise_args(args, bambu_arg_spec())) args$ndr <- 1 - testthat::expect_silent(bambu_validate_args(args)) + testthat::expect_silent(workflow_glue_r_normalise_args(args, bambu_arg_spec())) + + args$ndr <- NULL + args$threads <- "2" + normalised <- workflow_glue_r_normalise_args(args, bambu_arg_spec()) + testthat::expect_identical(normalised$threads, 2L) }) # Fail fast if --bams is empty rather than passing empty input to bambu. @@ -80,7 +87,7 @@ testthat::test_that("empty BAM list rejected", { ) testthat::expect_error( - bambu_resolve_inputs(args, bamfile_list_ctor = function(paths, yieldSize) paths), + bambu_resolve_inputs(args), "No BAM files were provided in --bams" ) }) @@ -94,13 +101,13 @@ testthat::test_that("unique sample aliases required", { sample_sheet = NULL ) testthat::expect_error( - bambu_resolve_inputs(args, bamfile_list_ctor = function(paths, yieldSize) paths), + bambu_resolve_inputs(args), "BAM aliases must be unique" ) args$aliases <- "sampleA" testthat::expect_error( - bambu_resolve_inputs(args, bamfile_list_ctor = function(paths, yieldSize) paths), + bambu_resolve_inputs(args), "Provide one alias per BAM in --bams" ) missing_alias_sheet <- tempfile(fileext = ".csv") @@ -118,7 +125,7 @@ testthat::test_that("unique sample aliases required", { sample_sheet = missing_alias_sheet ) testthat::expect_error( - bambu_resolve_inputs(args, bamfile_list_ctor = function(paths, yieldSize) paths), + bambu_resolve_inputs(args), "Sample sheet must contain an 'alias' column" ) @@ -134,7 +141,7 @@ testthat::test_that("unique sample aliases required", { ) args$sample_sheet <- duplicate_alias_sheet testthat::expect_error( - bambu_resolve_inputs(args, bamfile_list_ctor = function(paths, yieldSize) paths), + bambu_resolve_inputs(args), "Sample sheet aliases must be unique" ) }) @@ -159,13 +166,11 @@ testthat::test_that("sample sheet reordered to match BAMs", { aliases = "sampleA,sampleB", sample_sheet = sample_sheet ) - resolved <- bambu_resolve_inputs( - args, - bamfile_list_ctor = function(paths, yieldSize) paths - ) + resolved <- bambu_resolve_inputs(args) testthat::expect_equal(resolved$aliases, c("sampleA", "sampleB")) testthat::expect_equal(resolved$sample_df$alias, c("sampleA", "sampleB")) + testthat::expect_s4_class(resolved$reads, "BamFileList") bad_sheet <- tempfile(fileext = ".csv") writeLines( @@ -179,7 +184,7 @@ testthat::test_that("sample sheet reordered to match BAMs", { args$sample_sheet <- bad_sheet testthat::expect_error( - bambu_resolve_inputs(args, bamfile_list_ctor = function(paths, yieldSize) paths), + bambu_resolve_inputs(args), "Sample sheet is missing alias rows" ) }) @@ -199,6 +204,19 @@ testthat::test_that("transcriptome mode mapped to bambu args", { testthat::expect_true(discover$discovery) testthat::expect_equal(discover$NDR, 0.2) testthat::expect_equal(discover$ncore, 3L) + testthat::expect_true(discover$lowMemory) + testthat::expect_equal(discover$yieldSize, 250000L) + + auto_ndr_args <- list( + genome = "genome.fa", + threads = 2, + transcriptome_mode = "discover", + ndr = NULL + ) + auto_ndr <- bambu_build_args(auto_ndr_args, reads = "sample.bam", annotation_obj = annotation_obj) + testthat::expect_true(auto_ndr$discovery) + testthat::expect_false("NDR" %in% names(auto_ndr)) + testthat::expect_equal(auto_ndr$yieldSize, 250000L) fixed_args <- list( genome = "genome.fa", @@ -209,6 +227,8 @@ testthat::test_that("transcriptome mode mapped to bambu args", { fixed <- bambu_build_args(fixed_args, reads = "sample.bam", annotation_obj = annotation_obj) testthat::expect_false(fixed$discovery) testthat::expect_false("NDR" %in% names(fixed)) + testthat::expect_true(fixed$lowMemory) + testthat::expect_equal(fixed$yieldSize, 250000L) }) # End-to-end unit test with mocked bambu analysis function. @@ -235,19 +255,32 @@ testthat::test_that("bams input with discovery mode", { ) captured <- new.env(parent = emptyenv()) - fake_bamfile_list <- function(paths, yieldSize) { - captured$bamfile_paths <- paths - captured$yield_size <- yieldSize - structure(paths, class = "mockBamFileList") - } - fake_analysis <- function(reads, annotations, genome, ncore, discovery, NDR) { + fake_analysis <- function( + reads, + annotations, + genome, + ncore, + discovery, + lowMemory, + NDR = NULL, + yieldSize = NULL, + ... + ) { captured$reads <- reads captured$annotations <- annotations captured$genome <- genome captured$ncore <- ncore captured$discovery <- discovery captured$NDR <- NDR - make_test_tx_se(sample_names = c("sampleA", "sampleB")) + captured$low_memory <- lowMemory + captured$arg_yield_size <- yieldSize + captured$yield_size <- Rsamtools::yieldSize(reads[[1]]) + + base_se <- make_test_tx_se(sample_names = c("sampleA", "sampleB")) + SummarizedExperiment::SummarizedExperiment( + assays = SummarizedExperiment::assays(base_se), + rowRanges = make_test_bambu_row_ranges(fixture_dir) + ) } argv <- list( @@ -262,6 +295,7 @@ testthat::test_that("bams input with discovery mode", { threads = 2 ) + argv <- workflow_glue_r_normalise_args(argv, bambu_arg_spec()) result <- suppressMessages(main_run_bambu( argv, analysis_fn = fake_analysis, @@ -269,23 +303,17 @@ testthat::test_that("bams input with discovery mode", { captured$annotation_path <- annotation structure(list(path = annotation), class = "mockAnnotation") }, - gene_expression_fn = function(se) make_test_gene_se(sample_names = colnames(se)), - write_gtf_fn = function(row_ranges, file) { - writeLines( - 'chr1\tsim\texon\t1\t50\t.\t+\t.\tgene_id "gene1"; transcript_id "tx1";', - file - ) - }, - bamfile_list_ctor = fake_bamfile_list + gene_expression_fn = function(se) make_test_gene_se(sample_names = colnames(se)) )) testthat::expect_equal(captured$annotation_path, "annotation.gtf") testthat::expect_equal(captured$genome, "genome.fa") - testthat::expect_equal(captured$ncore, 2L) + testthat::expect_equal(captured$ncore, 1L) testthat::expect_true(captured$discovery) testthat::expect_equal(captured$NDR, 0.25) - testthat::expect_equal(captured$yield_size, 1000000) - testthat::expect_equal(captured$bamfile_paths, c(sample_a, sample_b)) + testthat::expect_true(captured$low_memory) + testthat::expect_equal(captured$arg_yield_size, 250000L) + testthat::expect_equal(captured$yield_size, 250000) testthat::expect_equal(result$sample_df$alias, c("sampleA", "sampleB")) testthat::expect_equal(result$sample_df$condition, c("control", "treated")) testthat::expect_true(file.exists(file.path(argv$out_dir, "bambu_qc_stats.json"))) @@ -308,35 +336,19 @@ testthat::test_that("list columns flattened for TSV output", { testthat::expect_equal(normalised$list_col, c("x;y", "z")) }) -# NCBI annotations may have gene_id="transcript_id" (literal string, not value reference). -# Replace with gene_id=, but leave gene_id="MYTRANSCRIPT_ID" alone. -testthat::test_that("malformed gene_id values sanitized", { - gtf_path <- tempfile(fileext = ".gtf") - writeLines( - c( - 'chr1\tsim\texon\t1\t50\t.\t+\t.\tgene_id "MYTRANSCRIPT_ID"; transcript_id "tx_keep";', - 'chr1\tsim\texon\t101\t150\t.\t+\t.\tgene_id "transcript_id"; transcript_id "tx_replace";' - ), - gtf_path - ) - - testthat::expect_warning( - bambu_sanitise_gtf_file(gtf_path), - "Replaced malformed gene_id" - ) - - lines <- readLines(gtf_path, warn = FALSE) - testthat::expect_match(lines[[1]], 'gene_id "MYTRANSCRIPT_ID";', fixed = TRUE) - testthat::expect_match(lines[[2]], 'gene_id "tx_replace";', fixed = TRUE) -}) - # Verify all expected output files are created with correct structure. testthat::test_that("bambu outputs written correctly", { out_dir <- tempfile("bambu-write-") dir.create(out_dir) sample_names <- c("sampleA", "sampleB") - se <- make_test_tx_se(sample_names = sample_names) + base_se <- make_test_tx_se(sample_names = sample_names) + row_ranges <- make_test_bambu_row_ranges(out_dir) + + se <- SummarizedExperiment::SummarizedExperiment( + assays = SummarizedExperiment::assays(base_se), + rowRanges = row_ranges + ) gene_se <- make_test_gene_se(sample_names = sample_names) sample_df <- data.frame(alias = sample_names, stringsAsFactors = FALSE) argv <- list(out_dir = out_dir, transcriptome_mode = "discover", ndr = 0.15) @@ -358,13 +370,7 @@ testthat::test_that("bambu outputs written correctly", { gene_se, sample_df, argv, - qc_stats, - write_gtf_fn = function(row_ranges, file) { - writeLines( - 'chr1\tsim\texon\t1\t50\t.\t+\t.\tgene_id "gene1"; transcript_id "tx1";', - file - ) - } + qc_stats ) testthat::expect_true(file.exists(file.path(out_dir, "transcripts.gtf"))) 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 263e487..ded16d2 100644 --- a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R +++ b/bin/workflow_glue_r/tests/testthat/test_de_analysis.R @@ -8,7 +8,7 @@ # Covariates CLI arg is comma-separated string with possible whitespace/empty values. testthat::test_that("covariates parsed and trimmed", { testthat::expect_equal( - de_parse_covariates(" batch, sex ,, site "), + workflow_glue_r_parse_csv_list(" batch, sex ,, site "), c("batch", "sex", "site") ) }) @@ -26,6 +26,7 @@ testthat::test_that("sample sheet structure validated", { covariates = "batch", reference_level = NULL ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) missing_alias <- data.frame(condition = rep(c("control", "treated"), each = 3)) testthat::expect_error( @@ -78,6 +79,7 @@ testthat::test_that("formula-unsafe column names rejected", { covariates = "batch", reference_level = NULL ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) testthat::expect_error( de_validate_inputs(tx_se, gene_se, sample_df, argv), @@ -150,6 +152,7 @@ testthat::test_that("non-syntactic aliases allowed, sheets reordered", { covariates = "batch", reference_level = "control" ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) validated <- de_validate_inputs(tx_se, gene_se, sample_df, argv) @@ -173,6 +176,7 @@ testthat::test_that("reference level defaults to control", { covariates = "batch", reference_level = NULL ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) sample_df <- data.frame( alias = colnames(tx_se), condition = rep(c("control", "treated"), each = 3), @@ -220,6 +224,7 @@ testthat::test_that("explicit reference level required without control", { covariates = "batch", reference_level = NULL ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) testthat::expect_error( de_validate_inputs(tx_se, gene_se, sample_df, argv), @@ -246,6 +251,7 @@ testthat::test_that("unusable count data rejected", { covariates = "batch", reference_level = "control" ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) sample_df <- data.frame( alias = colnames(tx_se), condition = rep(c("control", "treated"), each = 3), @@ -326,6 +332,7 @@ testthat::test_that("transcript SE without GENEID rejected", { reference_level = "control", out_dir = file.path(fixture_dir, "out") ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) testthat::expect_error( main_run_de_analysis(argv), @@ -369,6 +376,7 @@ testthat::test_that("underspecified designs rejected", { reference_level = "control", out_dir = file.path(fixture_dir, "out") ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) testthat::expect_error( suppressWarnings(main_run_de_analysis(argv)), @@ -621,6 +629,7 @@ testthat::test_that("contrasts expanded and samples subsetted", { out_dir = file.path(fixture_dir, "out") ) ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) suppressWarnings(suppressMessages(main_run_de_analysis(argv))) @@ -680,6 +689,7 @@ testthat::test_that("pipe characters in transcript IDs preserved", { out_dir = file.path(fixture_dir, "out") ) ) + argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) suppressWarnings(suppressMessages(main_run_de_analysis(argv))) diff --git a/subworkflows/transcriptome.nf b/subworkflows/transcriptome.nf index 4b9d380..f5eefd4 100644 --- a/subworkflows/transcriptome.nf +++ b/subworkflows/transcriptome.nf @@ -61,7 +61,7 @@ process runJointBambu { --transcriptome_mode "${params.transcriptome_mode}" \ --threads ${task.cpus} \ ${ndr_arg} \ - --out_dir cohort + --out_dir cohort \ """ }