kfukuta/wf-transcriptomes-v202
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

[CW-7214] Tidy up CLI parsing

Browse Source
This commit is contained in:
Chris Wright 2026-05-14 17:41:21 +00:00
parent a08e23849d
commit 5af948dabf
10 changed files with 674 additions and 443 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
.gitlab-ci.yml
View File

@ -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;"

9
bin/workflow_glue_r/NAMESPACE
View File

@ -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)

213
bin/workflow_glue_r/R/args.R Normal file
View File

@ -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
}

431
bin/workflow_glue_r/R/bambu.R
View File

@ -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 = ", ")
bambu_arg_spec <- function() {
list(
list(
name = "bams",
flag = "--bams",
help = "Comma-separated BAM paths.",
type = "character",
required = TRUE
),
call. = FALSE
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)
}

39
bin/workflow_glue_r/R/common.R
View File

@ -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_)

178
bin/workflow_glue_r/R/de_analysis.R
View File

@ -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",
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 <- 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) {
@ -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) {
" - <contrast>/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)
}

65
bin/workflow_glue_r/tests/testthat/helper-fixtures.R
View File

@ -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.

142
bin/workflow_glue_r/tests/testthat/test_bambu.R
View File

@ -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=<actual transcript_id value>, 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")))

12
bin/workflow_glue_r/tests/testthat/test_de_analysis.R
View File

@ -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)))

2
subworkflows/transcriptome.nf
View File

@ -61,7 +61,7 @@ process runJointBambu {
--transcriptome_mode "${params.transcriptome_mode}" \
--threads ${task.cpus} \
${ndr_arg} \
--out_dir cohort
--out_dir cohort \
"""
}