diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 8550b90..eb48389 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -13,6 +13,8 @@ variables: CI_FLAVOUR: "new" PYTEST_CONTAINER_NAME: "wf-common" PYTEST_CONTAINER_CONFIG_KEY: "common_sha" + RTEST_CONTAINER_NAME: "wf-transcriptomes-core" + RTEST_CONTAINER_CONFIG_KEY: "container_sha" macos-run: # Let's avoid those ARM64 runners for now diff --git a/bin/run_bambu.R b/bin/run_bambu.R deleted file mode 100755 index ab428ef..0000000 --- a/bin/run_bambu.R +++ /dev/null @@ -1,449 +0,0 @@ -#!/usr/bin/env Rscript - -# Set seed for reproducibility -set.seed(42) - -suppressPackageStartupMessages({ - library(argparser) - library(bambu) - library(Rsamtools) - library(SummarizedExperiment) - library(jsonlite) -}) - -parser <- arg_parser("Run bambu transcript discovery and quantification.") -parser <- add_argument(parser, "--bam_dir", help = "Directory containing BAM files.") -parser <- add_argument(parser, "--bam_path", help = "Path to a single BAM file.") -parser <- add_argument(parser, "--sample_alias", help = "Alias to use for a single BAM file.") -parser <- add_argument(parser, "--sample_sheet", help = "Optional sample sheet CSV.") -parser <- add_argument(parser, "--annotation", help = "Reference annotation GTF/GFF.") -parser <- add_argument(parser, "--genome", help = "Reference genome FASTA.") -parser <- add_argument(parser, "--transcriptome_mode", help = "discover or fixed_annotation.", default = "discover") -parser <- add_argument(parser, "--threads", help = "Number of worker threads.", type = "numeric", default = 1) -parser <- add_argument(parser, "--ndr", help = "Optional novel discovery rate.", type = "numeric") -parser <- add_argument(parser, "--out_dir", help = "Output directory.") -argv <- parse_args(parser) - -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 -} - -required_args <- c("annotation", "genome", "out_dir") -missing_args <- required_args[vapply(required_args, function(arg_name) { - value <- argv[[arg_name]] - arg_missing(value) -}, logical(1))] -if (length(missing_args) > 0) { - stop(sprintf( - "Missing required arguments: %s", - paste(sprintf("--%s", missing_args), collapse = ", ") - )) -} - -if (arg_missing(argv$bam_dir) == arg_missing(argv$bam_path)) { - stop("Provide exactly one of --bam_dir or --bam_path.") -} - -dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE) - -sample_df <- NULL -if (!arg_missing(argv$sample_sheet)) { - sample_df <- read.csv(argv$sample_sheet, check.names = FALSE, stringsAsFactors = FALSE) - if (!"alias" %in% names(sample_df)) { - stop("Sample sheet must contain an 'alias' column.") - } -} - -strip_alias <- function(path) { - name <- basename(path) - name <- sub("\\.aligned\\.sorted\\.bam$", "", name) - name <- tools::file_path_sans_ext(name) - name -} - -if (!arg_missing(argv$bam_dir)) { - bam_paths <- sort(list.files(argv$bam_dir, pattern = "\\.bam$", full.names = TRUE)) - if (length(bam_paths) < 1) { - stop("No BAM files were found in bam_dir.") - } - aliases <- vapply(bam_paths, strip_alias, character(1)) -} else { - bam_paths <- argv$bam_path - aliases <- if (!arg_missing(argv$sample_alias)) argv$sample_alias else strip_alias(argv$bam_path) -} - -if (!is.null(sample_df)) { - missing_aliases <- setdiff(aliases, sample_df$alias) - if (length(missing_aliases) > 0) { - stop(sprintf( - "Sample sheet is missing alias rows for BAM files: %s", - paste(missing_aliases, collapse = ", ") - )) - } - sample_df <- sample_df[match(aliases, sample_df$alias), , drop = FALSE] -} else { - sample_df <- data.frame(alias = aliases, stringsAsFactors = FALSE) -} - -annotation_obj <- prepareAnnotations(argv$annotation) -reads <- if (length(bam_paths) == 1) bam_paths else BamFileList(bam_paths, yieldSize = 1000000) - -# Handle NDR parameter with validation and documentation -default_ndr <- 0.1 -ndr_value <- default_ndr - -if (!arg_missing(argv$ndr)) { - if (argv$ndr < 0 || argv$ndr > 1) { - stop("NDR (Novel Discovery Rate) must be between 0 and 1") - } - ndr_value <- argv$ndr - message(sprintf("Using user-specified NDR = %.3f", ndr_value)) -} else { - message(sprintf("Using default NDR = %.3f", default_ndr)) -} - -if (identical(argv$transcriptome_mode, "discover")) { - 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)) -} - -bambu_args <- list( - reads = reads, - annotations = annotation_obj, - genome = argv$genome, - ncore = as.integer(argv$threads), - discovery = identical(argv$transcriptome_mode, "discover") -) - -if (identical(argv$transcriptome_mode, "discover")) { - bambu_args$NDR <- ndr_value -} - -message("Running bambu...") -se <- do.call(bambu, bambu_args) -message("Bambu completed successfully") -colnames(se) <- aliases - -counts_mat <- assays(se)$counts -full_length_mat <- assays(se)$fullLengthCounts - -# Collect QC statistics before filtering -qc_stats <- list() -qc_stats$total_transcripts_before_filter <- nrow(se) -qc_stats$total_genes_before_filter <- length(unique(rowData(se)$GENEID)) -qc_stats$samples <- ncol(se) - -# Filter low-count transcripts -if (is.null(full_length_mat)) { - keep_idx <- rowSums(counts_mat) > 0 -} else { - keep_idx <- rowSums(full_length_mat) > 0 -} -if (!any(keep_idx)) { - keep_idx <- rowSums(counts_mat) >= 0 -} - -qc_stats$transcripts_filtered <- sum(!keep_idx) -message(sprintf("Filtering: keeping %d / %d transcripts", sum(keep_idx), length(keep_idx))) - -se <- se[keep_idx, ] - -# Library size statistics and warnings -lib_sizes <- colSums(assays(se)$counts) -qc_stats$library_sizes <- as.list(lib_sizes) -qc_stats$min_library_size <- min(lib_sizes) -qc_stats$max_library_size <- max(lib_sizes) -qc_stats$median_library_size <- median(lib_sizes) - -if (length(lib_sizes) > 1) { - lib_size_ratio <- max(lib_sizes) / min(lib_sizes) - qc_stats$library_size_ratio <- lib_size_ratio - - if (lib_size_ratio > 3) { - warning(sprintf( - "Large library size variation detected (%.1fx difference).\n Min: %d, Max: %d reads.\n CPM normalization may not be appropriate for such variation.", - lib_size_ratio, min(lib_sizes), max(lib_sizes) - )) - qc_stats$library_size_warning <- sprintf("%.1fx variation (>3x threshold)", lib_size_ratio) - } -} - -# Per-sample detection statistics -qc_stats$transcripts_detected_per_sample <- as.list(colSums(assays(se)$counts > 0)) -qc_stats$median_transcripts_detected <- median(colSums(assays(se)$counts > 0)) - -qc_stats$total_transcripts_after_filter <- nrow(se) -qc_stats$total_genes_after_filter <- length(unique(rowData(se)$GENEID)) - -row_ranges <- rowRanges(se) -writeToGTF(row_ranges, file = file.path(argv$out_dir, "transcripts.gtf")) - -gene_se <- transcriptToGeneExpression(se) -colnames(gene_se) <- aliases - -saveRDS(se, file.path(argv$out_dir, "bambu_transcripts.rds")) -saveRDS(gene_se, file.path(argv$out_dir, "bambu_genes.rds")) -write.csv(sample_df, file.path(argv$out_dir, "samples.csv"), row.names = FALSE, quote = FALSE) - -tx_meta <- as.data.frame(rowData(se)) -if (!"TXNAME" %in% names(tx_meta)) { - tx_meta$TXNAME <- rownames(se) -} -if (!"GENEID" %in% names(tx_meta)) { - tx_meta$GENEID <- NA_character_ -} -gene_meta <- as.data.frame(rowData(gene_se)) -if (!"GENEID" %in% names(gene_meta)) { - gene_meta$GENEID <- rownames(gene_se) -} - -matrix_to_df <- function(se_obj, assay_name, id_col, meta_df) { - assay_df <- as.data.frame(assays(se_obj)[[assay_name]]) - assay_df[[id_col]] <- rownames(se_obj) - assay_df <- assay_df[, c(id_col, setdiff(names(assay_df), id_col)), drop = FALSE] - merge(meta_df, assay_df, by.x = id_col, by.y = id_col, all.y = TRUE, sort = FALSE) -} - -normalise_tsv_value <- function(value) { - if (length(value) == 0 || all(is.na(value))) { - return(NA_character_) - } - if (is.list(value)) { - value <- unlist(value, recursive = TRUE, use.names = FALSE) - } - if (length(value) == 0 || all(is.na(value))) { - return(NA_character_) - } - paste(as.character(value), collapse = ";") -} - -normalise_tsv_df <- function(df) { - as.data.frame( - lapply(df, function(column) { - if (is.list(column)) { - vapply(column, normalise_tsv_value, character(1)) - } else { - column - } - }), - stringsAsFactors = FALSE, - check.names = FALSE - ) -} - -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] -} - -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 -} - -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 <- normalise_gtf_attribute_value( - extract_gtf_attribute(attr_field, "transcript_id") - ) - gene_id <- normalise_gtf_attribute_value( - extract_gtf_attribute(attr_field, "gene_id") - ) - - if (!is.null(gene_id) && grepl("\\btranscript_id\\b", gene_id)) { - 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) -} - -tx_meta <- normalise_tsv_df(tx_meta) -gene_meta <- normalise_tsv_df(gene_meta) - -write.table( - tx_meta, - file = file.path(argv$out_dir, "transcript_metadata.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE -) -write.table( - gene_meta, - file = file.path(argv$out_dir, "gene_metadata.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE -) - -tx_counts <- matrix_to_df(se, "counts", "TXNAME", tx_meta) -tx_cpm <- matrix_to_df(se, "CPM", "TXNAME", tx_meta) -gene_counts <- matrix_to_df(gene_se, "counts", "GENEID", gene_meta) -gene_cpm <- matrix_to_df(gene_se, "CPM", "GENEID", gene_meta) - -write.table( - tx_counts, - file = file.path(argv$out_dir, "transcript_counts.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE -) -write.table( - tx_cpm, - file = file.path(argv$out_dir, "transcript_cpm.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE -) -write.table( - gene_counts, - file = file.path(argv$out_dir, "gene_counts.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE -) -write.table( - gene_cpm, - file = file.path(argv$out_dir, "gene_cpm.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE -) - -sanitise_gtf_file(file.path(argv$out_dir, "transcripts.gtf")) - -# Write QC statistics as JSON for HTML report -qc_stats$transcriptome_mode <- argv$transcriptome_mode -qc_stats$ndr_used <- if (identical(argv$transcriptome_mode, "discover")) ndr_value else "N/A" -qc_stats$timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S") - -write_json( - qc_stats, - file.path(argv$out_dir, "bambu_qc_stats.json"), - pretty = TRUE, - auto_unbox = TRUE -) - -format_count <- function(value) { - if (length(value) == 0 || all(is.na(value))) { - return("NA") - } - format( - round(as.numeric(value), 0), - scientific = FALSE, - trim = TRUE, - big.mark = "," - ) -} - -# Write human-readable QC summary -qc_summary <- c( - "Bambu Quantification QC Summary", - "================================", - "", - sprintf("Timestamp: %s", qc_stats$timestamp), - sprintf("Mode: %s", argv$transcriptome_mode), - if (identical(argv$transcriptome_mode, "discover")) sprintf("NDR: %.3f", ndr_value) else NULL, - "", - "Sample Statistics:", - sprintf(" Samples analyzed: %s", format_count(qc_stats$samples)), - sprintf( - " Median library size: %s reads", - format_count(qc_stats$median_library_size) - ), - sprintf( - " Library size range: %s - %s reads", - format_count(qc_stats$min_library_size), - format_count(qc_stats$max_library_size) - ), - if (!is.null(qc_stats$library_size_warning)) sprintf(" WARNING: %s", qc_stats$library_size_warning) else NULL, - "", - "Transcript Discovery:", - sprintf( - " Transcripts before filtering: %s", - format_count(qc_stats$total_transcripts_before_filter) - ), - sprintf( - " Transcripts after filtering: %s", - format_count(qc_stats$total_transcripts_after_filter) - ), - sprintf( - " Transcripts removed: %s", - format_count(qc_stats$transcripts_filtered) - ), - sprintf( - " Median transcripts detected per sample: %s", - format_count(qc_stats$median_transcripts_detected) - ), - "", - "Gene-Level Summary:", - sprintf( - " Unique genes (before filter): %s", - format_count(qc_stats$total_genes_before_filter) - ), - sprintf( - " Unique genes (after filter): %s", - format_count(qc_stats$total_genes_after_filter) - ), - "" -) - -writeLines(qc_summary, file.path(argv$out_dir, "bambu_qc_summary.txt")) -message("QC statistics written to bambu_qc_stats.json and bambu_qc_summary.txt") - -# Save session info for reproducibility -writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt")) -message("Session info saved for reproducibility") diff --git a/bin/run_de_analysis.R b/bin/run_de_analysis.R deleted file mode 100755 index 9bdd12b..0000000 --- a/bin/run_de_analysis.R +++ /dev/null @@ -1,717 +0,0 @@ -#!/usr/bin/env Rscript - -# Set seed for reproducibility -set.seed(42) - -suppressPackageStartupMessages({ - library(argparser) - library(DESeq2) - library(DEXSeq) - library(SummarizedExperiment) - library(jsonlite) -}) - -parser <- arg_parser("Run DESeq2 and DEXSeq on bambu output.") -parser <- add_argument(parser, "--transcript_rds", help = "bambu transcript RDS.") -parser <- add_argument(parser, "--gene_rds", help = "bambu gene RDS.") -parser <- add_argument(parser, "--sample_sheet", help = "Sample sheet CSV.") -parser <- add_argument(parser, "--condition_column", help = "Primary condition column.", default = "condition") -parser <- add_argument(parser, "--covariates", help = "Comma-separated nuisance covariates.") -parser <- add_argument(parser, "--reference_level", help = "Reference level for the condition column.") -parser <- add_argument(parser, "--out_dir", help = "Output directory.", default = "de_analysis") -argv <- parse_args(parser) - -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 -} - -required_args <- c("transcript_rds", "gene_rds", "sample_sheet") -missing_args <- required_args[vapply(required_args, function(arg_name) { - value <- argv[[arg_name]] - arg_missing(value) -}, logical(1))] -if (length(missing_args) > 0) { - stop(sprintf( - "Missing required arguments: %s", - paste(sprintf("--%s", missing_args), collapse = ", ") - )) -} - -dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE) - -tx_se <- readRDS(argv$transcript_rds) -gene_se <- readRDS(argv$gene_rds) -sample_df <- read.csv(argv$sample_sheet, check.names = FALSE, stringsAsFactors = FALSE) - -if (!"alias" %in% names(sample_df)) { - stop("Sample sheet must contain an 'alias' column.") -} -if (!(argv$condition_column %in% names(sample_df))) { - stop(sprintf("Sample sheet must contain the '%s' column.", argv$condition_column)) -} - -covariates <- character(0) -if (!arg_missing(argv$covariates)) { - covariates <- trimws(strsplit(argv$covariates, ",", fixed = TRUE)[[1]]) - covariates <- covariates[nzchar(covariates)] -} -missing_covariates <- setdiff(covariates, names(sample_df)) -if (length(missing_covariates) > 0) { - stop(sprintf("Missing covariate columns: %s", paste(missing_covariates, collapse = ", "))) -} - -sample_df <- sample_df[match(colnames(tx_se), sample_df$alias), , drop = FALSE] -if (any(is.na(sample_df$alias))) { - stop("Sample sheet aliases do not match the bambu output sample names.") -} - -condition_values <- unique(sample_df[[argv$condition_column]]) -if (length(condition_values) < 2) { - stop("Differential analysis requires at least two condition levels.") -} - -reference_level <- argv$reference_level -if (arg_missing(reference_level)) { - if ("control" %in% condition_values) { - reference_level <- "control" - } else { - stop("Provide --reference_level when the condition column does not contain 'control'.") - } -} -if (!(reference_level %in% condition_values)) { - stop("The requested reference level is not present in the condition column.") -} - -sample_df[[argv$condition_column]] <- factor(sample_df[[argv$condition_column]]) -for (covariate in covariates) { - sample_df[[covariate]] <- factor(sample_df[[covariate]]) -} - -run_deseq_with_fallback <- function(dds, contrast_name = "unknown") { - tryCatch( - DESeq(dds, quiet = TRUE), - error = function(err) { - if (!grepl( - "all gene-wise dispersion estimates are within 2 orders of magnitude", - conditionMessage(err), - fixed = TRUE - )) { - stop(err) - } - - warning( - "STATISTICAL POWER REDUCED: DESeq2 dispersion estimation failed for ", contrast_name, ".\n", - "This usually indicates:\n", - " 1. Too few replicates (recommend n>=3 per group)\n", - " 2. High biological variability\n", - " 3. Poor data quality\n", - "Falling back to gene-wise dispersion (no information sharing).\n", - "Results will have reduced power and wider confidence intervals." - ) - - dds <- estimateSizeFactors(dds) - dds <- estimateDispersionsGeneEst(dds) - dispersions(dds) <- mcols(dds)$dispGeneEst - - # Write diagnostic file - diag_content <- c( - "DESeq2 Dispersion Estimation Fallback Applied", - "==============================================", - "", - sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), - sprintf("Contrast: %s", contrast_name), - sprintf("Samples: %d", ncol(dds)), - sprintf("Genes tested: %d", nrow(dds)), - sprintf("Dispersion range: %.3f to %.3f", min(dispersions(dds)), max(dispersions(dds))), - "", - "WHAT HAPPENED:", - " Curve fitting failed. Using gene-wise dispersion estimates.", - "", - "IMPLICATIONS:", - " - No information sharing across genes", - " - Reduced statistical power", - " - Wider confidence intervals", - " - More conservative results (fewer discoveries)", - "", - "LIKELY CAUSES:", - " 1. Too few replicates (recommend n>=3 per group)", - " 2. High biological variability", - " 3. Poor data quality or outlier samples", - "", - "RECOMMENDATIONS:", - " - Add more biological replicates if possible", - " - Check sample quality metrics", - " - Consider filtering low-count genes more stringently" - ) - - diag_file <- file.path(argv$out_dir, sprintf("DESeq2_dispersion_fallback_%s.txt", gsub("[^A-Za-z0-9_-]", "_", contrast_name))) - writeLines(diag_content, diag_file) - - nbinomWaldTest(dds) - } - ) -} - -estimate_dispersions_with_fallback <- function(object, context_label, allow_gene_est = TRUE) { - tryCatch( - estimateDispersions(object), - error = function(err) { - if (!grepl( - "all gene-wise dispersion estimates are within 2 orders of magnitude", - conditionMessage(err), - fixed = TRUE - )) { - stop(err) - } - - message( - context_label, - " dispersion fitting failed; ", - "retrying with fitType='local'." - ) - tryCatch( - estimateDispersions(object, fitType = "local"), - error = function(local_err) { - if (!grepl( - "all gene-wise dispersion estimates are within 2 orders of magnitude", - conditionMessage(local_err), - fixed = TRUE - )) { - stop(local_err) - } - - message( - context_label, - " local-fit dispersion retry failed; ", - "retrying with fitType='mean'." - ) - tryCatch( - estimateDispersions(object, fitType = "mean"), - error = function(mean_err) { - if (!grepl( - "all gene-wise dispersion estimates are within 2 orders of magnitude", - conditionMessage(mean_err), - fixed = TRUE - )) { - stop(mean_err) - } - if (!allow_gene_est) { - stop(mean_err) - } - - message( - context_label, - " mean-fit dispersion retry failed; ", - "falling back to gene-wise dispersion estimates." - ) - object <- estimateDispersionsGeneEst(object) - dispersions(object) <- mcols(object)$dispGeneEst - object - } - ) - } - ) - } - ) -} - -normalise_tsv_value <- function(value) { - if (length(value) == 0 || all(is.na(value))) { - return(NA_character_) - } - if (is.list(value)) { - value <- unlist(value, recursive = TRUE, use.names = FALSE) - } - if (length(value) == 0 || all(is.na(value))) { - return(NA_character_) - } - paste(as.character(value), collapse = ";") -} - -normalise_tsv_df <- function(df) { - as.data.frame( - lapply(df, function(column) { - if (is.list(column)) { - vapply(column, normalise_tsv_value, character(1)) - } else { - column - } - }), - stringsAsFactors = FALSE, - check.names = FALSE - ) -} - -is_recoverable_dexseq_error <- function(message_text) { - grepl( - "all gene-wise dispersion estimates are within 2 orders of magnitude", - message_text, - fixed = TRUE - ) || grepl( - "model matrix is not full rank", - message_text, - fixed = TRUE - ) || grepl( - "replacement has 1 row, data has 0", - message_text, - fixed = TRUE - ) -} - -empty_tsv <- function(columns) { - out <- as.data.frame(matrix(nrow = 0, ncol = length(columns))) - names(out) <- columns - out -} - -write_placeholder_pdf <- function(path, label) { - pdf(path) - plot.new() - text(0.5, 0.5, label, cex = 0.9) - dev.off() -} - -run_deseq2 <- function(count_mat, coldata, target_level, contrast_name) { - design_terms <- c(covariates, argv$condition_column) - design_formula <- as.formula(paste("~", paste(design_terms, collapse = " + "))) - dds <- DESeqDataSetFromMatrix( - countData = round(count_mat), - colData = coldata, - design = design_formula - ) - dds <- run_deseq_with_fallback(dds, contrast_name) - results(dds, contrast = c(argv$condition_column, target_level, reference_level), independentFiltering = TRUE) -} - -run_dexseq <- function(tx_counts, tx_meta, coldata, active_covariates = covariates) { - coldata$sample <- factor(coldata$alias) - coldata[[argv$condition_column]] <- factor(coldata[[argv$condition_column]]) - for (covariate in active_covariates) { - coldata[[covariate]] <- factor(coldata[[covariate]]) - } - - covariate_exon_terms <- if (length(active_covariates) > 0) { - paste0(active_covariates, ":exon") - } else { - character(0) - } - design_terms <- c("sample", "exon", covariate_exon_terms, paste0(argv$condition_column, ":exon")) - reduced_terms <- c("sample", "exon", covariate_exon_terms) - full_formula <- as.formula(paste("~", paste(design_terms, collapse = " + "))) - reduced_formula <- as.formula(paste("~", paste(reduced_terms, collapse = " + "))) - - tryCatch({ - dxd <- DEXSeqDataSet( - countData = round(tx_counts), - sampleData = as.data.frame(coldata), - design = full_formula, - featureID = tx_meta$TXNAME, - groupID = tx_meta$GENEID - ) - dxd <- estimateSizeFactors(dxd) - dxd <- estimate_dispersions_with_fallback(dxd, "DEXSeq", allow_gene_est = TRUE) - dxd <- testForDEU(dxd, reducedModel = reduced_formula) - dxd <- estimateExonFoldChanges(dxd, fitExpToVar = argv$condition_column) - dxr <- DEXSeqResults(dxd, independentFiltering = FALSE) - list(dxd = dxd, dxr = dxr) - }, error = function(err) { - if (length(active_covariates) == 0 || !grepl( - "model matrix is not full rank", - conditionMessage(err), - fixed = TRUE - )) { - stop(err) - } - - dropped_covariate <- tail(active_covariates, 1) - kept_covariates <- head(active_covariates, -1) - message( - "DEXSeq design was not full rank with covariate '", - dropped_covariate, - "'; retrying without it." - ) - run_dexseq(tx_counts, tx_meta, coldata, kept_covariates) - }) -} - -tx_meta <- as.data.frame(rowData(tx_se)) -if (!"TXNAME" %in% names(tx_meta)) { - tx_meta$TXNAME <- rownames(tx_se) -} -if (!"GENEID" %in% names(tx_meta)) { - stop("Transcript rowData must contain GENEID for DEXSeq.") -} - -gene_meta <- as.data.frame(rowData(gene_se)) -if (!"GENEID" %in% names(gene_meta)) { - gene_meta$GENEID <- rownames(gene_se) -} - -targets <- setdiff(as.character(condition_values), reference_level) - -# Initialize QC statistics collector -de_qc_stats <- list() -de_qc_stats$timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S") -de_qc_stats$total_samples <- nrow(sample_df) -de_qc_stats$condition_column <- argv$condition_column -de_qc_stats$reference_level <- reference_level -de_qc_stats$covariates <- if (length(covariates) > 0) covariates else "none" -de_qc_stats$num_contrasts <- length(targets) -de_qc_stats$contrasts <- list() - -# Check sample sizes and warn if underpowered -n_per_group <- table(sample_df[[argv$condition_column]]) -de_qc_stats$samples_per_group <- as.list(n_per_group) - -sample_size_warnings <- c() -if (any(n_per_group < 3)) { - warning( - "WARNING: Some condition groups have fewer than 3 replicates.\n", - "Recommended minimum for DGE: n=3 per group\n", - "Current sample sizes: ", paste(names(n_per_group), "=", n_per_group, collapse=", "), "\n", - "Results may have reduced statistical power." - ) - sample_size_warnings <- c(sample_size_warnings, "Some groups have n<3 (recommended minimum)") -} - -if (any(n_per_group < 2)) { - stop("ERROR: Some condition groups have fewer than 2 replicates. Cannot perform statistical testing.") -} - -de_qc_stats$sample_size_warnings <- if (length(sample_size_warnings) > 0) sample_size_warnings else "none" - -# Multiple testing warning -if (length(targets) > 1) { - fwer <- (1 - (1-0.05)^length(targets)) * 100 - mt_warning <- sprintf( - "Multiple contrasts tested (%d). Per-contrast FDR < 0.05 yields family-wise error rate of ~%.1f%%", - length(targets), fwer - ) - message("WARNING: ", mt_warning) - de_qc_stats$multiple_testing_note <- mt_warning - - mt_content <- c( - "Multiple Testing Across Contrasts", - "==================================", - "", - sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), - sprintf("Number of contrasts tested: %d", length(targets)), - sprintf("Contrasts: %s", paste(sprintf("%s vs %s", targets, reference_level), collapse=", ")), - "", - "PER-CONTRAST FDR THRESHOLD: 0.05", - sprintf("FAMILY-WISE ERROR RATE: ~%.1f%%", fwer), - "", - "WHAT THIS MEANS:", - " Each contrast uses FDR < 0.05 independently.", - " When testing multiple contrasts, the overall false positive rate increases.", - sprintf(" Expected: %.1f%% chance of at least one false positive across all contrasts", fwer), - "", - "RECOMMENDATIONS:", - " 1. Use stricter per-contrast threshold:", - sprintf(" Bonferroni correction: 0.05 / %d = %.4f", length(targets), 0.05/length(targets)), - " 2. Focus on pre-specified contrasts of interest", - " 3. Treat results as exploratory and validate key findings", - " 4. Consider using hierarchical testing procedures", - "", - "INTERPRETATION:", - " - Results passing FDR < 0.05 in each contrast are discoveries for that contrast", - " - But the overall false discovery burden is higher than 5%", - " - Prioritize genes significant across multiple contrasts", - " - Validate top findings experimentally" - ) - writeLines(mt_content, file.path(argv$out_dir, "MULTIPLE_TESTING_WARNING.txt")) -} - -for (target_level in targets) { - contrast_name <- sprintf("%s_%s_vs_%s", argv$condition_column, target_level, reference_level) - contrast_dir <- file.path(argv$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) - contrast_samples <- droplevels(sample_df[keep_samples, , drop = FALSE]) - contrast_samples[[argv$condition_column]] <- relevel( - factor(contrast_samples[[argv$condition_column]]), - ref = reference_level - ) - - # Collect per-contrast QC stats - contrast_qc <- list() - contrast_qc$name <- contrast_name - contrast_qc$target_level <- target_level - contrast_qc$reference_level <- reference_level - contrast_qc$n_samples <- nrow(contrast_samples) - contrast_qc$n_target <- sum(contrast_samples[[argv$condition_column]] == target_level) - contrast_qc$n_reference <- sum(contrast_samples[[argv$condition_column]] == reference_level) - - # DTU power warning - if (nrow(contrast_samples) < 6) { - dtu_warning <- sprintf( - "DTU analysis may be underpowered (n=%d, recommend n>=6 with >=3 per group)", - nrow(contrast_samples) - ) - warning(dtu_warning) - contrast_qc$dtu_power_warning <- dtu_warning - } - - gene_counts <- assays(gene_se)$counts[, contrast_samples$alias, drop = FALSE] - tx_counts <- assays(tx_se)$counts[, contrast_samples$alias, drop = FALSE] - - contrast_qc$genes_tested <- nrow(gene_counts) - contrast_qc$transcripts_tested <- nrow(tx_counts) - - dge_res <- as.data.frame(run_deseq2(gene_counts, contrast_samples, target_level, contrast_name)) - dge_res$GENEID <- rownames(dge_res) - dge_res <- merge(gene_meta, dge_res, by = "GENEID", all.y = TRUE, sort = FALSE) - dge_res <- normalise_tsv_df(dge_res) - - # Collect DGE statistics - contrast_qc$dge_total_genes <- nrow(dge_res) - contrast_qc$dge_significant_fdr05 <- sum(dge_res$padj < 0.05, na.rm = TRUE) - contrast_qc$dge_significant_fdr01 <- sum(dge_res$padj < 0.01, na.rm = TRUE) - contrast_qc$dge_upregulated <- sum(dge_res$padj < 0.05 & dge_res$log2FoldChange > 0, na.rm = TRUE) - contrast_qc$dge_downregulated <- sum(dge_res$padj < 0.05 & dge_res$log2FoldChange < 0, na.rm = TRUE) - - write.table( - dge_res[order(dge_res$padj), ], - file = file.path(contrast_dir, "results_dge.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE - ) - - pdf(file.path(contrast_dir, "results_dge.pdf")) - dds_plot <- DESeqDataSetFromMatrix( - countData = round(gene_counts), - colData = contrast_samples, - design = as.formula(paste("~", paste(c(covariates, argv$condition_column), collapse = " + "))) - ) - dds_plot <- run_deseq_with_fallback(dds_plot, contrast_name) - plotMA(results(dds_plot, contrast = c(argv$condition_column, target_level, reference_level), independentFiltering = TRUE)) - dev.off() - - dex_res <- tryCatch( - run_dexseq(tx_counts, tx_meta, contrast_samples), - error = function(err) { - message_text <- conditionMessage(err) - if (!is_recoverable_dexseq_error(message_text)) { - stop(err) - } - - warning( - "DEXSeq failed for contrast ", target_level, " vs ", reference_level, "\n", - "Error: ", message_text - ) - - # Write explicit failure report - failure_content <- c( - "DTU Analysis Failed", - "===================", - "", - sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), - sprintf("Contrast: %s vs %s", target_level, reference_level), - sprintf("Samples: %d (%d %s, %d %s)", - nrow(contrast_samples), - sum(contrast_samples[[argv$condition_column]] == target_level), target_level, - sum(contrast_samples[[argv$condition_column]] == reference_level), reference_level), - sprintf("Transcripts: %d", nrow(tx_counts)), - "", - "ERROR MESSAGE:", - sprintf(" %s", message_text), - "", - "DTU RESULTS CANNOT BE INTERPRETED", - "", - "This failure is likely due to:", - " 1. Insufficient samples (need >=3 per group, recommend >=6 total for DTU)", - " 2. Too few transcripts with sufficient counts", - " 3. Design matrix not full rank (covariate confounding)", - " 4. Extreme count distributions", - "", - "RECOMMENDATIONS:", - " - Use gene-level DGE results (less power required)", - " - Add more biological replicates", - " - Filter transcripts more stringently", - " - Simplify experimental design (remove problematic covariates)", - "", - "NOTE: Empty DTU result files indicate analysis failure, not 'no DTU detected'" - ) - writeLines(failure_content, file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt")) - - NULL - } - ) - - if (is.null(dex_res)) { - dex_df <- empty_tsv(c( - "featureID", - "groupID", - "log2fold", - "pvalue", - "padj", - "exonBaseMean" - )) - tx_dtu <- dex_df - gene_dtu <- empty_tsv(c("GENEID", "qval")) - write_placeholder_pdf( - file.path(contrast_dir, "results_dtu.pdf"), - "DEXSeq did not converge for this contrast.\nSee DTU_ANALYSIS_FAILED.txt for details." - ) - contrast_qc$dtu_status <- "FAILED" - contrast_qc$dtu_significant_transcripts <- 0 - contrast_qc$dtu_significant_genes <- 0 - } else { - dxr <- dex_res$dxr - dxd <- dex_res$dxd - dex_df <- as.data.frame(dxr) - dex_df <- normalise_tsv_df(dex_df) - tx_dtu <- dex_df[, intersect( - c("featureID", "groupID", "log2fold", "pvalue", "padj", "exonBaseMean"), - names(dex_df) - ), drop = FALSE] - tx_dtu <- normalise_tsv_df(tx_dtu) - - gene_q <- perGeneQValue(dxr) - gene_dtu <- data.frame( - GENEID = names(gene_q), - qval = unname(gene_q), - row.names = NULL - ) - - # Collect DTU statistics - contrast_qc$dtu_status <- "SUCCESS" - contrast_qc$dtu_significant_transcripts <- sum(tx_dtu$padj < 0.05, na.rm = TRUE) - contrast_qc$dtu_significant_genes <- sum(gene_dtu$qval < 0.05, na.rm = TRUE) - - pdf(file.path(contrast_dir, "results_dtu.pdf")) - plotMA(dxr, cex = 0.8, alpha = 0.05) - plotDispEsts(dxd) - dev.off() - } - - write.table( - dex_df, - file = file.path(contrast_dir, "results_dexseq.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE - ) - - write.table( - tx_dtu[order(tx_dtu$padj), ], - file = file.path(contrast_dir, "results_dtu_transcript.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE - ) - - write.table( - gene_dtu[order(gene_dtu$qval), ], - file = file.path(contrast_dir, "results_dtu_gene.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE - ) - - write.table( - contrast_samples, - file = file.path(contrast_dir, "samples_used.tsv"), - sep = "\t", - quote = FALSE, - row.names = FALSE - ) - - # Write per-contrast QC summary - contrast_qc_summary <- c( - sprintf("Contrast QC Summary: %s", contrast_name), - paste(rep("=", 50), collapse = ""), - "", - "Sample Information:", - sprintf(" Target level (%s): %d samples", target_level, contrast_qc$n_target), - sprintf(" Reference level (%s): %d samples", reference_level, contrast_qc$n_reference), - sprintf(" Total samples: %d", contrast_qc$n_samples), - "", - "DGE Results:", - sprintf(" Genes tested: %d", contrast_qc$genes_tested), - sprintf(" Significant (FDR < 0.05): %d", contrast_qc$dge_significant_fdr05), - sprintf(" Significant (FDR < 0.01): %d", contrast_qc$dge_significant_fdr01), - sprintf(" Upregulated: %d", contrast_qc$dge_upregulated), - sprintf(" Downregulated: %d", contrast_qc$dge_downregulated), - "", - "DTU Results:", - sprintf(" Status: %s", contrast_qc$dtu_status), - sprintf(" Transcripts tested: %d", contrast_qc$transcripts_tested), - if (contrast_qc$dtu_status == "SUCCESS") { - c( - sprintf(" Significant transcripts (FDR < 0.05): %d", contrast_qc$dtu_significant_transcripts), - sprintf(" Genes with DTU (q < 0.05): %d", contrast_qc$dtu_significant_genes) - ) - } else { - " See DTU_ANALYSIS_FAILED.txt for details" - }, - if (!is.null(contrast_qc$dtu_power_warning)) paste0(" WARNING: ", contrast_qc$dtu_power_warning) else NULL, - "" - ) - writeLines(contrast_qc_summary, file.path(contrast_dir, "contrast_qc_summary.txt")) - - # Add to overall QC stats - de_qc_stats$contrasts[[contrast_name]] <- contrast_qc -} - -# Write overall DE/DTU QC statistics as JSON for HTML report -write_json( - de_qc_stats, - file.path(argv$out_dir, "de_qc_stats.json"), - pretty = TRUE, - auto_unbox = TRUE -) - -# Write human-readable overall summary -overall_summary <- c( - "Differential Expression/Usage Analysis Summary", - paste(rep("=", 50), collapse = ""), - "", - sprintf("Timestamp: %s", de_qc_stats$timestamp), - sprintf("Total samples: %d", de_qc_stats$total_samples), - sprintf("Condition column: %s", de_qc_stats$condition_column), - sprintf("Reference level: %s", de_qc_stats$reference_level), - sprintf("Covariates: %s", paste(de_qc_stats$covariates, collapse = ", ")), - "", - "Sample Sizes:", - sapply(names(de_qc_stats$samples_per_group), function(grp) { - sprintf(" %s: %d samples", grp, de_qc_stats$samples_per_group[[grp]]) - }), - if (de_qc_stats$sample_size_warnings != "none") paste0(" WARNING: ", de_qc_stats$sample_size_warnings) else NULL, - "", - sprintf("Number of contrasts tested: %d", de_qc_stats$num_contrasts), - if (!is.null(de_qc_stats$multiple_testing_note)) paste0(" NOTE: ", de_qc_stats$multiple_testing_note) else NULL, - "", - "Per-Contrast Results:", - sapply(names(de_qc_stats$contrasts), function(cname) { - cqc <- de_qc_stats$contrasts[[cname]] - c( - "", - sprintf(" %s:", cname), - sprintf(" Samples: %d (%d vs %d)", cqc$n_samples, cqc$n_target, cqc$n_reference), - sprintf(" DGE significant: %d genes (FDR<0.05)", cqc$dge_significant_fdr05), - sprintf(" DTU status: %s", cqc$dtu_status), - if (cqc$dtu_status == "SUCCESS") sprintf(" DTU significant: %d genes", cqc$dtu_significant_genes) else NULL - ) - }), - "", - "For detailed per-contrast statistics, see:", - " - /contrast_qc_summary.txt", - " - /results_dge.tsv", - " - /results_dtu_gene.tsv", - "" -) -writeLines(unlist(overall_summary), file.path(argv$out_dir, "de_overall_summary.txt")) - -# Save session info for reproducibility -writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt")) -message("QC statistics written to de_qc_stats.json and de_overall_summary.txt") -message("Session info saved for reproducibility") diff --git a/bin/supeRglue b/bin/supeRglue new file mode 100755 index 0000000..e4377ec --- /dev/null +++ b/bin/supeRglue @@ -0,0 +1,24 @@ +#!/usr/bin/env Rscript + +supeRglue_script_dir <- function() { + file_arg <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE) + if (length(file_arg) > 0) { + return(dirname(normalizePath(sub("^--file=", "", file_arg[[1]]), mustWork = TRUE))) + } + + frame_files <- vapply(sys.frames(), function(frame) { + path <- frame$ofile + if (is.null(path)) "" else path + }, character(1)) + frame_files <- frame_files[nzchar(frame_files)] + if (length(frame_files) > 0) { + return(dirname(normalizePath(frame_files[[length(frame_files)]], mustWork = TRUE))) + } + + stop("Unable to locate supeRglue on disk.", call. = FALSE) +} + +script_dir <- supeRglue_script_dir() +source(file.path(script_dir, "workflow_glue_r", "load.R")) +workflow_glue_r_load(file.path(script_dir, "workflow_glue_r")) +workflow_glue_r_cli() diff --git a/bin/workflow_glue_r/DESCRIPTION b/bin/workflow_glue_r/DESCRIPTION new file mode 100644 index 0000000..c4ce4f3 --- /dev/null +++ b/bin/workflow_glue_r/DESCRIPTION @@ -0,0 +1,20 @@ +Package: workflowGlueR +Title: Local R Helpers for wf-transcriptomes +Version: 0.0.1.0000 +Authors@R: + person("Oxford Nanopore Technologies", role = c("aut", "cre")) +Description: Local R helper functions used by wf-transcriptomes. +License: Proprietary +Encoding: UTF-8 +LazyData: false +Imports: + argparser, + bambu, + DESeq2, + DEXSeq, + jsonlite, + Rsamtools, + S4Vectors, + SummarizedExperiment +Suggests: + testthat diff --git a/bin/workflow_glue_r/NAMESPACE b/bin/workflow_glue_r/NAMESPACE new file mode 100644 index 0000000..6d935f5 --- /dev/null +++ b/bin/workflow_glue_r/NAMESPACE @@ -0,0 +1,26 @@ +export(bambu_arg_parser) +export(bambu_build_args) +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_cli) +export(workflow_glue_r_components) +export(workflow_glue_r_arg_missing) +export(workflow_glue_r_empty_tsv) +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/bambu.R b/bin/workflow_glue_r/R/bambu.R new file mode 100644 index 0000000..1512201 --- /dev/null +++ b/bin/workflow_glue_r/R/bambu.R @@ -0,0 +1,538 @@ +bambu_arg_parser <- function() { + parser <- argparser::arg_parser("Run bambu transcript discovery and quantification.") + parser <- argparser::add_argument(parser, "--bam_dir", help = "Directory containing BAM files.") + parser <- argparser::add_argument(parser, "--bam_path", help = "Path to a single BAM file.") + parser <- argparser::add_argument(parser, "--sample_alias", help = "Alias to use for a single BAM file.") + 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$bam_dir) == workflow_glue_r_arg_missing(argv$bam_path)) { + stop("Provide exactly one of --bam_dir or --bam_path.", 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 + ) + } + + 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_strip_alias <- function(path) { + name <- basename(path) + name <- sub("\\.aligned\\.sorted\\.bam$", "", name) + tools::file_path_sans_ext(name) +} + +bambu_resolve_inputs <- function( + argv, + bamfile_list_ctor = Rsamtools::BamFileList, + list_files_fn = base::list.files +) { + sample_df <- NULL + if (!workflow_glue_r_arg_missing(argv$sample_sheet)) { + sample_df <- workflow_glue_r_read_csv(argv$sample_sheet) + if (!"alias" %in% names(sample_df)) { + stop("Sample sheet must contain an 'alias' column.", call. = FALSE) + } + duplicate_sample_aliases <- unique(sample_df$alias[duplicated(sample_df$alias)]) + if (length(duplicate_sample_aliases) > 0) { + stop( + sprintf( + "Sample sheet aliases must be unique; duplicated aliases: %s", + paste(duplicate_sample_aliases, collapse = ", ") + ), + call. = FALSE + ) + } + } + + if (!workflow_glue_r_arg_missing(argv$bam_dir)) { + bam_paths <- sort(list_files_fn(argv$bam_dir, pattern = "\\.bam$", full.names = TRUE)) + if (length(bam_paths) < 1) { + stop("No BAM files were found in bam_dir.", call. = FALSE) + } + aliases <- unname(vapply(bam_paths, bambu_strip_alias, character(1))) + } else { + bam_paths <- argv$bam_path + aliases <- if (!workflow_glue_r_arg_missing(argv$sample_alias)) { + argv$sample_alias + } else { + bambu_strip_alias(argv$bam_path) + } + } + + duplicate_bam_aliases <- unique(aliases[duplicated(aliases)]) + if (length(duplicate_bam_aliases) > 0) { + stop( + sprintf( + "BAM aliases must be unique; duplicated aliases: %s", + paste(duplicate_bam_aliases, collapse = ", ") + ), + call. = FALSE + ) + } + + if (!is.null(sample_df)) { + missing_aliases <- setdiff(aliases, sample_df$alias) + if (length(missing_aliases) > 0) { + stop( + sprintf( + "Sample sheet is missing alias rows for BAM files: %s", + paste(missing_aliases, collapse = ", ") + ), + call. = FALSE + ) + } + sample_df <- sample_df[match(aliases, sample_df$alias), , drop = FALSE] + } else { + sample_df <- data.frame(alias = aliases, stringsAsFactors = FALSE) + } + + reads <- if (length(bam_paths) == 1) { + bam_paths + } else { + bamfile_list_ctor(bam_paths, yieldSize = 1000000) + } + + list( + bam_paths = bam_paths, + aliases = aliases, + sample_df = sample_df, + reads = reads + ) +} + +bambu_discovery_enabled <- function(argv) { + identical(argv$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_args <- list( + reads = reads, + annotations = annotation_obj, + genome = argv$genome, + ncore = as.integer(argv$threads), + discovery = bambu_discovery_enabled(argv) + ) + + if (bambu_discovery_enabled(argv)) { + bambu_args$NDR <- bambu_resolve_ndr(argv) + } + + bambu_args +} + +bambu_filter_transcripts <- function(se) { + counts_mat <- SummarizedExperiment::assays(se)$counts + full_length_mat <- SummarizedExperiment::assays(se)$fullLengthCounts + + gene_ids <- SummarizedExperiment::rowData(se)$GENEID + qc_stats <- list( + total_transcripts_before_filter = nrow(se), + total_genes_before_filter = length(unique(gene_ids)), + samples = ncol(se) + ) + + if (is.null(full_length_mat)) { + keep_idx <- rowSums(counts_mat) > 0 + } else { + keep_idx <- rowSums(full_length_mat) > 0 + } + if (!any(keep_idx)) { + keep_idx <- rowSums(counts_mat) >= 0 + } + + qc_stats$transcripts_filtered <- sum(!keep_idx) + se <- se[keep_idx, ] + qc_stats$total_transcripts_after_filter <- nrow(se) + qc_stats$total_genes_after_filter <- length(unique(SummarizedExperiment::rowData(se)$GENEID)) + + list(se = se, qc_stats = qc_stats) +} + +bambu_matrix_to_df <- function(se_obj, assay_name, id_col, meta_df) { + assay_df <- as.data.frame(SummarizedExperiment::assays(se_obj)[[assay_name]]) + assay_df[[id_col]] <- rownames(se_obj) + assay_df <- assay_df[, c(id_col, setdiff(names(assay_df), id_col)), drop = FALSE] + 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") + } + format( + round(as.numeric(value), 0), + scientific = FALSE, + trim = TRUE, + big.mark = "," + ) +} + +bambu_write_outputs <- function(se, gene_se, sample_df, argv, qc_stats, write_gtf_fn = bambu::writeToGTF) { + write_gtf_fn( + SummarizedExperiment::rowRanges(se), + file = file.path(argv$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")) + utils::write.csv( + sample_df, + file.path(argv$out_dir, "samples.csv"), + row.names = FALSE, + quote = FALSE + ) + + tx_meta <- as.data.frame(SummarizedExperiment::rowData(se)) + if (!"TXNAME" %in% names(tx_meta)) { + tx_meta$TXNAME <- rownames(se) + } + if (!"GENEID" %in% names(tx_meta)) { + tx_meta$GENEID <- NA_character_ + } + + gene_meta <- as.data.frame(SummarizedExperiment::rowData(gene_se)) + if (!"GENEID" %in% names(gene_meta)) { + gene_meta$GENEID <- rownames(gene_se) + } + + tx_meta <- workflow_glue_r_normalise_tsv_df(tx_meta) + gene_meta <- workflow_glue_r_normalise_tsv_df(gene_meta) + + utils::write.table( + tx_meta, + file = file.path(argv$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"), + 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 + ) + utils::write.table( + tx_cpm, + file = file.path(argv$out_dir, "transcript_cpm.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + utils::write.table( + gene_counts, + file = file.path(argv$out_dir, "gene_counts.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + utils::write.table( + gene_cpm, + file = file.path(argv$out_dir, "gene_cpm.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + + 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) + } else { + "N/A" + } + qc_stats$timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S") + + jsonlite::write_json( + qc_stats, + file.path(argv$out_dir, "bambu_qc_stats.json"), + pretty = TRUE, + auto_unbox = TRUE + ) + + qc_summary <- c( + "Bambu Quantification QC Summary", + "================================", + "", + 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, + "", + "Sample Statistics:", + sprintf(" Samples analyzed: %s", bambu_format_count(qc_stats$samples)), + sprintf(" Median library size: %s reads", bambu_format_count(qc_stats$median_library_size)), + sprintf( + " Library size range: %s - %s reads", + bambu_format_count(qc_stats$min_library_size), + bambu_format_count(qc_stats$max_library_size) + ), + if (!is.null(qc_stats$library_size_warning)) sprintf(" WARNING: %s", qc_stats$library_size_warning) else NULL, + "", + "Transcript Discovery:", + sprintf(" Transcripts before filtering: %s", bambu_format_count(qc_stats$total_transcripts_before_filter)), + sprintf(" Transcripts after filtering: %s", bambu_format_count(qc_stats$total_transcripts_after_filter)), + sprintf(" Transcripts removed: %s", bambu_format_count(qc_stats$transcripts_filtered)), + sprintf( + " Median transcripts detected per sample: %s", + bambu_format_count(qc_stats$median_transcripts_detected) + ), + "", + "Gene-Level Summary:", + sprintf(" Unique genes (before filter): %s", bambu_format_count(qc_stats$total_genes_before_filter)), + sprintf(" Unique genes (after filter): %s", bambu_format_count(qc_stats$total_genes_after_filter)), + "" + ) + + writeLines(qc_summary, file.path(argv$out_dir, "bambu_qc_summary.txt")) + writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt")) +} + +main_run_bambu <- function( + argv, + analysis_fn = bambu::bambu, + prepare_annotations_fn = bambu::prepareAnnotations, + gene_expression_fn = bambu::transcriptToGeneExpression, + write_gtf_fn = bambu::writeToGTF, + bamfile_list_ctor = Rsamtools::BamFileList, + list_files_fn = base::list.files +) { + set.seed(42) + suppressPackageStartupMessages({ + library(GenomicRanges) + library(Rsamtools) + }) + + bambu_validate_args(argv) + dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE) + + inputs <- bambu_resolve_inputs( + argv, + bamfile_list_ctor = bamfile_list_ctor, + list_files_fn = list_files_fn + ) + annotation_obj <- prepare_annotations_fn(argv$annotation) + ndr_value <- bambu_resolve_ndr(argv) + + 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 (bambu_discovery_enabled(argv)) { + 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)) + } + + message("Running bambu...") + se <- do.call(analysis_fn, bambu_build_args(argv, inputs$reads, annotation_obj)) + message("Bambu completed successfully") + colnames(se) <- inputs$aliases + + filtered <- bambu_filter_transcripts(se) + se <- filtered$se + qc_stats <- filtered$qc_stats + message( + sprintf( + "Filtering: keeping %d / %d transcripts", + qc_stats$total_transcripts_after_filter, + qc_stats$total_transcripts_before_filter + ) + ) + + lib_sizes <- colSums(SummarizedExperiment::assays(se)$counts) + qc_stats$library_sizes <- as.list(lib_sizes) + qc_stats$min_library_size <- min(lib_sizes) + qc_stats$max_library_size <- max(lib_sizes) + qc_stats$median_library_size <- stats::median(lib_sizes) + + if (length(lib_sizes) > 1) { + lib_size_ratio <- max(lib_sizes) / min(lib_sizes) + qc_stats$library_size_ratio <- lib_size_ratio + if (lib_size_ratio > 3) { + warning( + sprintf( + paste0( + "Large library size variation detected (%.1fx difference).\n", + " Min: %d, Max: %d reads.\n", + " CPM normalization may not be appropriate for such variation." + ), + lib_size_ratio, + min(lib_sizes), + max(lib_sizes) + ) + ) + qc_stats$library_size_warning <- sprintf("%.1fx variation (>3x threshold)", lib_size_ratio) + } + } + + detected_per_sample <- colSums(SummarizedExperiment::assays(se)$counts > 0) + qc_stats$transcripts_detected_per_sample <- as.list(detected_per_sample) + qc_stats$median_transcripts_detected <- stats::median(detected_per_sample) + + gene_se <- gene_expression_fn(se) + colnames(gene_se) <- inputs$aliases + + bambu_write_outputs( + se, + gene_se, + inputs$sample_df, + argv, + qc_stats, + write_gtf_fn = write_gtf_fn + ) + + invisible( + list( + se = se, + gene_se = gene_se, + sample_df = inputs$sample_df, + qc_stats = qc_stats + ) + ) +} + +run_bambu_cli <- function(argv = commandArgs(trailingOnly = TRUE)) { + parsed <- argparser::parse_args(bambu_arg_parser(), argv = argv) + main_run_bambu(parsed) +} diff --git a/bin/workflow_glue_r/R/cli.R b/bin/workflow_glue_r/R/cli.R new file mode 100644 index 0000000..abc8de4 --- /dev/null +++ b/bin/workflow_glue_r/R/cli.R @@ -0,0 +1,72 @@ +workflow_glue_r_components <- function(env = globalenv()) { + parser_suffix <- "_arg_parser" + parser_names <- grep( + paste0(parser_suffix, "$"), + ls(envir = env, all.names = TRUE), + value = TRUE + ) + + components <- list() + for (parser_name in parser_names) { + component <- sub(paste0(parser_suffix, "$"), "", parser_name) + cli_name <- paste0("run_", component, "_cli") + parser <- get(parser_name, envir = env) + runner <- if (exists(cli_name, envir = env, mode = "function")) { + get(cli_name, envir = env) + } else { + NULL + } + + if (is.function(parser) && is.function(runner)) { + components[[component]] <- list( + name = component, + parser_name = parser_name, + runner_name = cli_name, + parser = parser, + runner = runner + ) + } + } + + components[sort(names(components))] +} + +workflow_glue_r_usage <- function(components = workflow_glue_r_components()) { + component_names <- names(components) + lines <- c( + "Usage: supeRglue [options]", + "", + "Commands:", + if (length(component_names) > 0) { + paste0(" ", component_names) + } else { + " " + }, + "", + "Use 'supeRglue --help' for command-specific options." + ) + paste(lines, collapse = "\n") +} + +workflow_glue_r_cli <- function(argv = commandArgs(trailingOnly = TRUE), env = globalenv()) { + components <- workflow_glue_r_components(env = env) + + if (length(argv) < 1 || argv[[1]] %in% c("-h", "--help", "help")) { + cat(workflow_glue_r_usage(components), "\n") + return(invisible(0L)) + } + + command <- argv[[1]] + if (!command %in% names(components)) { + stop( + sprintf( + "Unknown supeRglue command '%s'. Available commands: %s", + command, + paste(names(components), collapse = ", ") + ), + call. = FALSE + ) + } + + components[[command]]$runner(argv[-1]) +} diff --git a/bin/workflow_glue_r/R/common.R b/bin/workflow_glue_r/R/common.R new file mode 100644 index 0000000..0fa0caa --- /dev/null +++ b/bin/workflow_glue_r/R/common.R @@ -0,0 +1,94 @@ +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 && + grepl("^[A-Za-z][A-Za-z0-9_.]*$", name) +} + +workflow_glue_r_validate_r_formula_names <- function(names, label = "Column") { + invalid <- names[!vapply(names, workflow_glue_r_is_r_formula_name, logical(1))] + if (length(invalid) > 0) { + stop( + sprintf( + "%s names must be safe for R formulas. Invalid names: %s. Names must start with a letter and contain only letters, numbers, underscores, and dots.", + label, + paste(invalid, collapse = ", ") + ), + call. = FALSE + ) + } + 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_) + } + if (is.list(value)) { + value <- unlist(value, recursive = TRUE, use.names = FALSE) + } + if (length(value) == 0 || all(is.na(value))) { + return(NA_character_) + } + paste(as.character(value), collapse = ";") +} + +workflow_glue_r_normalise_tsv_df <- function(df) { + as.data.frame( + lapply(df, function(column) { + if (is.list(column)) { + vapply(column, workflow_glue_r_normalise_tsv_value, character(1)) + } else { + column + } + }), + stringsAsFactors = FALSE, + check.names = FALSE + ) +} + +bambu_normalise_tsv_df <- workflow_glue_r_normalise_tsv_df + +workflow_glue_r_empty_tsv <- function(columns) { + out <- as.data.frame(matrix(nrow = 0, ncol = length(columns))) + names(out) <- columns + out +} diff --git a/bin/workflow_glue_r/R/de_analysis.R b/bin/workflow_glue_r/R/de_analysis.R new file mode 100644 index 0000000..ee6a140 --- /dev/null +++ b/bin/workflow_glue_r/R/de_analysis.R @@ -0,0 +1,781 @@ +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" + ) + 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_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) + workflow_glue_r_validate_r_formula_names( + c(argv$condition_column, covariates), + label = "Design column" + ) + + if (!"alias" %in% names(sample_df)) { + stop("Sample sheet must contain an 'alias' column.", call. = FALSE) + } + duplicate_sample_aliases <- unique(sample_df$alias[duplicated(sample_df$alias)]) + if (length(duplicate_sample_aliases) > 0) { + stop( + sprintf( + "Sample sheet aliases must be unique; duplicated aliases: %s", + paste(duplicate_sample_aliases, collapse = ", ") + ), + call. = FALSE + ) + } + if (!(argv$condition_column %in% names(sample_df))) { + stop( + sprintf("Sample sheet must contain the '%s' column.", argv$condition_column), + call. = FALSE + ) + } + + missing_covariates <- setdiff(covariates, names(sample_df)) + if (length(missing_covariates) > 0) { + stop( + sprintf( + "Missing covariate columns: %s", + paste(missing_covariates, collapse = ", ") + ), + call. = FALSE + ) + } + + if (any(duplicated(colnames(tx_se)))) { + stop("Transcript RDS sample names must be unique.", call. = FALSE) + } + if (any(duplicated(colnames(gene_se)))) { + stop("Gene RDS sample names must be unique.", call. = FALSE) + } + if (!setequal(colnames(tx_se), colnames(gene_se))) { + stop("Transcript and gene RDS sample names must match.", call. = FALSE) + } + + tx_counts <- SummarizedExperiment::assays(tx_se)$counts + gene_counts <- SummarizedExperiment::assays(gene_se)$counts + if (is.null(tx_counts)) { + stop("Transcript RDS must contain a 'counts' assay.", call. = FALSE) + } + if (is.null(gene_counts)) { + stop("Gene RDS must contain a 'counts' assay.", call. = FALSE) + } + if (anyNA(tx_counts) || anyNA(gene_counts)) { + stop("Count matrices must not contain NA values.", call. = FALSE) + } + zero_count_samples <- unique(c( + colnames(tx_counts)[colSums(tx_counts) == 0], + colnames(gene_counts)[colSums(gene_counts) == 0] + )) + if (length(zero_count_samples) > 0) { + stop( + sprintf( + "Count matrices contain samples with zero total counts: %s", + paste(zero_count_samples, collapse = ", ") + ), + call. = FALSE + ) + } + + sample_df <- sample_df[match(colnames(tx_se), sample_df$alias), , drop = FALSE] + if (any(is.na(sample_df$alias))) { + stop("Sample sheet aliases do not match the bambu output sample names.", call. = FALSE) + } + + condition_values <- unique(sample_df[[argv$condition_column]]) + if (length(condition_values) < 2) { + stop("Differential analysis requires at least two condition levels.", call. = FALSE) + } + + reference_level <- argv$reference_level + if (workflow_glue_r_arg_missing(reference_level)) { + if ("control" %in% condition_values) { + reference_level <- "control" + } else { + stop( + "Provide --reference_level when the condition column does not contain 'control'.", + call. = FALSE + ) + } + } + if (!(reference_level %in% condition_values)) { + stop("The requested reference level is not present in the condition column.", call. = FALSE) + } + + sample_df[[argv$condition_column]] <- factor(sample_df[[argv$condition_column]]) + for (covariate in covariates) { + sample_df[[covariate]] <- factor(sample_df[[covariate]]) + } + + list( + tx_se = tx_se, + gene_se = gene_se, + sample_df = sample_df, + covariates = covariates, + condition_values = condition_values, + reference_level = reference_level + ) +} + +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_run_deseq_with_fallback <- function(dds, contrast_name, out_dir) { + tryCatch( + DESeq2::DESeq(dds, quiet = TRUE), + error = function(err) { + if (!grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + conditionMessage(err), + fixed = TRUE + )) { + stop(err) + } + + warning( + "STATISTICAL POWER REDUCED: DESeq2 dispersion estimation failed for ", + contrast_name, + ".\n", + "This usually indicates:\n", + " 1. Too few replicates (recommend n>=3 per group)\n", + " 2. High biological variability\n", + " 3. Poor data quality\n", + "Falling back to gene-wise dispersion (no information sharing).\n", + "Results will have reduced power and wider confidence intervals." + ) + + dds <- DESeq2::estimateSizeFactors(dds) + dds <- DESeq2::estimateDispersionsGeneEst(dds) + dds <- de_set_dispersions(dds, S4Vectors::mcols(dds)$dispGeneEst) + + diag_content <- c( + "DESeq2 Dispersion Estimation Fallback Applied", + "==============================================", + "", + sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), + sprintf("Contrast: %s", contrast_name), + sprintf("Samples: %d", ncol(dds)), + sprintf("Genes tested: %d", nrow(dds)), + sprintf( + "Dispersion range: %.3f to %.3f", + min(DESeq2::dispersions(dds)), + max(DESeq2::dispersions(dds)) + ), + "", + "WHAT HAPPENED:", + " Curve fitting failed. Using gene-wise dispersion estimates.", + "", + "IMPLICATIONS:", + " - No information sharing across genes", + " - Reduced statistical power", + " - Wider confidence intervals", + " - More conservative results (fewer discoveries)", + "", + "LIKELY CAUSES:", + " 1. Too few replicates (recommend n>=3 per group)", + " 2. High biological variability", + " 3. Poor data quality or outlier samples", + "", + "RECOMMENDATIONS:", + " - Add more biological replicates if possible", + " - Check sample quality metrics", + " - Consider filtering low-count genes more stringently" + ) + + diag_file <- file.path( + out_dir, + sprintf( + "DESeq2_dispersion_fallback_%s.txt", + gsub("[^A-Za-z0-9_-]", "_", contrast_name) + ) + ) + writeLines(diag_content, diag_file) + + DESeq2::nbinomWaldTest(dds) + } + ) +} + +de_estimate_dispersions_with_fallback <- function(object, context_label, allow_gene_est = TRUE) { + tryCatch( + DESeq2::estimateDispersions(object), + error = function(err) { + if (!grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + conditionMessage(err), + fixed = TRUE + )) { + stop(err) + } + + message(context_label, " dispersion fitting failed; retrying with fitType='local'.") + tryCatch( + DESeq2::estimateDispersions(object, fitType = "local"), + error = function(local_err) { + if (!grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + conditionMessage(local_err), + fixed = TRUE + )) { + stop(local_err) + } + + message(context_label, " local-fit dispersion retry failed; retrying with fitType='mean'.") + tryCatch( + DESeq2::estimateDispersions(object, fitType = "mean"), + error = function(mean_err) { + if (!grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + conditionMessage(mean_err), + fixed = TRUE + )) { + stop(mean_err) + } + if (!allow_gene_est) { + stop(mean_err) + } + + message( + context_label, + " mean-fit dispersion retry failed; falling back to gene-wise dispersion estimates." + ) + object <- DESeq2::estimateDispersionsGeneEst(object) + object <- de_set_dispersions(object, S4Vectors::mcols(object)$dispGeneEst) + object + } + ) + } + ) + } + ) +} + +de_is_recoverable_dexseq_error <- function(message_text) { + grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + message_text, + fixed = TRUE + ) || grepl( + "model matrix is not full rank", + message_text, + fixed = TRUE + ) || grepl( + "replacement has 1 row, data has 0", + message_text, + fixed = TRUE + ) +} + +de_write_placeholder_pdf <- function(path, label) { + grDevices::pdf(path) + graphics::plot.new() + graphics::text(0.5, 0.5, label, cex = 0.9) + grDevices::dev.off() +} + +de_run_deseq2_result <- function( + count_mat, + coldata, + target_level, + reference_level, + condition_column, + covariates, + out_dir, + contrast_name +) { + design_terms <- c(covariates, condition_column) + design_formula <- stats::as.formula(paste("~", paste(design_terms, collapse = " + "))) + dds <- DESeq2::DESeqDataSetFromMatrix( + countData = round(count_mat), + colData = coldata, + design = design_formula + ) + dds <- de_run_deseq_with_fallback(dds, contrast_name, out_dir) + result <- DESeq2::results( + dds, + contrast = c(condition_column, target_level, reference_level), + independentFiltering = TRUE + ) + list(dds = dds, result = result) +} + +de_run_dexseq_result <- function( + tx_counts, + tx_meta, + coldata, + condition_column, + covariates +) { + coldata$sample <- factor(coldata$alias) + coldata[[condition_column]] <- factor(coldata[[condition_column]]) + for (covariate in covariates) { + coldata[[covariate]] <- factor(coldata[[covariate]]) + } + + run_inner <- function(active_covariates) { + covariate_exon_terms <- if (length(active_covariates) > 0) { + paste0(active_covariates, ":exon") + } else { + character(0) + } + design_terms <- c("sample", "exon", covariate_exon_terms, paste0(condition_column, ":exon")) + reduced_terms <- c("sample", "exon", covariate_exon_terms) + full_formula <- stats::as.formula(paste("~", paste(design_terms, collapse = " + "))) + reduced_formula <- stats::as.formula(paste("~", paste(reduced_terms, collapse = " + "))) + + tryCatch({ + dxd <- DEXSeq::DEXSeqDataSet( + countData = round(tx_counts), + sampleData = as.data.frame(coldata), + design = full_formula, + featureID = tx_meta$TXNAME, + groupID = tx_meta$GENEID + ) + dxd <- DESeq2::estimateSizeFactors(dxd) + dxd <- de_estimate_dispersions_with_fallback(dxd, "DEXSeq", allow_gene_est = TRUE) + dxd <- DEXSeq::testForDEU(dxd, reducedModel = reduced_formula) + dxd <- DEXSeq::estimateExonFoldChanges(dxd, fitExpToVar = condition_column) + dxr <- DEXSeq::DEXSeqResults(dxd, independentFiltering = FALSE) + list(dxd = dxd, dxr = dxr) + }, error = function(err) { + if (length(active_covariates) == 0 || !grepl( + "model matrix is not full rank", + conditionMessage(err), + fixed = TRUE + )) { + stop(err) + } + + dropped_covariate <- tail(active_covariates, 1) + kept_covariates <- head(active_covariates, -1) + message( + "DEXSeq design was not full rank with covariate '", + dropped_covariate, + "'; retrying without it." + ) + run_inner(kept_covariates) + }) + } + + run_inner(covariates) +} + +main_run_de_analysis <- function( + argv, + deseq_runner = de_run_deseq2_result, + dexseq_runner = de_run_dexseq_result, + pdf_fn = grDevices::pdf, + dev_off_fn = grDevices::dev.off, + plot_ma_fn = DESeq2::plotMA, + plot_disp_fn = DESeq2::plotDispEsts, + per_gene_q_fn = DEXSeq::perGeneQValue, + placeholder_pdf_fn = de_write_placeholder_pdf +) { + set.seed(42) + dir.create(argv$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) + sample_df <- validated$sample_df + covariates <- validated$covariates + condition_values <- validated$condition_values + reference_level <- validated$reference_level + + tx_meta <- as.data.frame(SummarizedExperiment::rowData(tx_se)) + if (!"TXNAME" %in% names(tx_meta)) { + tx_meta$TXNAME <- rownames(tx_se) + } + if (!"GENEID" %in% names(tx_meta)) { + stop("Transcript rowData must contain GENEID for DEXSeq.", call. = FALSE) + } + + gene_meta <- as.data.frame(SummarizedExperiment::rowData(gene_se)) + if (!"GENEID" %in% names(gene_meta)) { + gene_meta$GENEID <- rownames(gene_se) + } + + targets <- setdiff(as.character(condition_values), reference_level) + + de_qc_stats <- list( + timestamp = format(Sys.time(), "%Y-%m-%d %H:%M:%S"), + total_samples = nrow(sample_df), + condition_column = argv$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]]) + de_qc_stats$samples_per_group <- as.list(n_per_group) + + sample_size_warnings <- character(0) + if (any(n_per_group < 3)) { + warning( + "WARNING: Some condition groups have fewer than 3 replicates.\n", + "Recommended minimum for DGE: n=3 per group\n", + "Current sample sizes: ", + paste(names(n_per_group), "=", n_per_group, collapse = ", "), + "\nResults may have reduced statistical power." + ) + sample_size_warnings <- c(sample_size_warnings, "Some groups have n<3 (recommended minimum)") + } + if (any(n_per_group < 2)) { + stop( + "ERROR: Some condition groups have fewer than 2 replicates. Cannot perform statistical testing.", + call. = FALSE + ) + } + de_qc_stats$sample_size_warnings <- if (length(sample_size_warnings) > 0) { + sample_size_warnings + } else { + "none" + } + + if (length(targets) > 1) { + fwer <- (1 - (1 - 0.05)^length(targets)) * 100 + mt_warning <- sprintf( + "Multiple contrasts tested (%d). Per-contrast FDR < 0.05 yields family-wise error rate of ~%.1f%%", + length(targets), + fwer + ) + message("WARNING: ", mt_warning) + de_qc_stats$multiple_testing_note <- mt_warning + mt_content <- c( + "Multiple Testing Across Contrasts", + "==================================", + "", + sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), + sprintf("Number of contrasts tested: %d", length(targets)), + sprintf("Contrasts: %s", paste(sprintf("%s vs %s", targets, reference_level), collapse = ", ")), + "", + "PER-CONTRAST FDR THRESHOLD: 0.05", + sprintf("FAMILY-WISE ERROR RATE: ~%.1f%%", fwer), + "", + "WHAT THIS MEANS:", + " Each contrast uses FDR < 0.05 independently.", + " When testing multiple contrasts, the overall false positive rate increases.", + sprintf(" Expected: %.1f%% chance of at least one false positive across all contrasts", fwer), + "", + "RECOMMENDATIONS:", + sprintf(" 1. Bonferroni correction: 0.05 / %d = %.4f", length(targets), 0.05 / length(targets)), + " 2. Focus on pre-specified contrasts of interest", + " 3. Treat results as exploratory and validate key findings", + " 4. Consider using hierarchical testing procedures", + "" + ) + writeLines(mt_content, file.path(argv$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) + dir.create(contrast_dir, showWarnings = FALSE, recursive = TRUE) + + keep_samples <- sample_df[[argv$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]]), + ref = reference_level + ) + + contrast_qc <- list( + name = contrast_name, + 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) + ) + + if (nrow(contrast_samples) < 6) { + contrast_qc$dtu_power_warning <- sprintf( + "DTU analysis may be underpowered (n=%d, recommend n>=6 with >=3 per group)", + nrow(contrast_samples) + ) + warning(contrast_qc$dtu_power_warning) + } + + gene_counts <- SummarizedExperiment::assays(gene_se)$counts[, contrast_samples$alias, drop = FALSE] + tx_counts <- SummarizedExperiment::assays(tx_se)$counts[, contrast_samples$alias, drop = FALSE] + contrast_qc$genes_tested <- nrow(gene_counts) + contrast_qc$transcripts_tested <- nrow(tx_counts) + + dge_run <- deseq_runner( + gene_counts, + contrast_samples, + target_level, + reference_level, + argv$condition_column, + covariates, + argv$out_dir, + contrast_name + ) + dge_res <- as.data.frame(dge_run$result) + dge_res$GENEID <- rownames(dge_res) + dge_res <- merge(gene_meta, dge_res, by = "GENEID", all.y = TRUE, sort = FALSE) + dge_res <- workflow_glue_r_normalise_tsv_df(dge_res) + + contrast_qc$dge_total_genes <- nrow(dge_res) + contrast_qc$dge_significant_fdr05 <- sum(dge_res$padj < 0.05, na.rm = TRUE) + contrast_qc$dge_significant_fdr01 <- sum(dge_res$padj < 0.01, na.rm = TRUE) + contrast_qc$dge_upregulated <- sum( + dge_res$padj < 0.05 & dge_res$log2FoldChange > 0, + na.rm = TRUE + ) + contrast_qc$dge_downregulated <- sum( + dge_res$padj < 0.05 & dge_res$log2FoldChange < 0, + na.rm = TRUE + ) + + utils::write.table( + dge_res[order(dge_res$padj), ], + file = file.path(contrast_dir, "results_dge.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + + pdf_fn(file.path(contrast_dir, "results_dge.pdf")) + plot_ma_fn(dge_run$result) + dev_off_fn() + + dex_res <- tryCatch( + dexseq_runner( + tx_counts, + tx_meta, + contrast_samples, + argv$condition_column, + covariates + ), + error = function(err) { + message_text <- conditionMessage(err) + if (!de_is_recoverable_dexseq_error(message_text)) { + stop(err) + } + + warning( + "DEXSeq failed for contrast ", + target_level, + " vs ", + reference_level, + "\nError: ", + message_text + ) + + failure_content <- c( + "DTU Analysis Failed", + "===================", + "", + sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), + sprintf("Contrast: %s vs %s", target_level, reference_level), + sprintf( + "Samples: %d (%d %s, %d %s)", + nrow(contrast_samples), + sum(contrast_samples[[argv$condition_column]] == target_level), + target_level, + sum(contrast_samples[[argv$condition_column]] == reference_level), + reference_level + ), + sprintf("Transcripts: %d", nrow(tx_counts)), + "", + "ERROR MESSAGE:", + sprintf(" %s", message_text), + "", + "DTU RESULTS CANNOT BE INTERPRETED", + "" + ) + writeLines(failure_content, file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt")) + NULL + } + ) + + if (is.null(dex_res)) { + dex_df <- workflow_glue_r_empty_tsv(c( + "featureID", + "groupID", + "log2fold", + "pvalue", + "padj", + "exonBaseMean" + )) + tx_dtu <- dex_df + gene_dtu <- workflow_glue_r_empty_tsv(c("GENEID", "qval")) + placeholder_pdf_fn( + file.path(contrast_dir, "results_dtu.pdf"), + "DEXSeq did not converge for this contrast.\nSee DTU_ANALYSIS_FAILED.txt for details." + ) + contrast_qc$dtu_status <- "FAILED" + contrast_qc$dtu_significant_transcripts <- 0 + contrast_qc$dtu_significant_genes <- 0 + } else { + dex_df <- as.data.frame(dex_res$dxr) + dex_df <- workflow_glue_r_normalise_tsv_df(dex_df) + tx_dtu <- dex_df[, intersect( + c("featureID", "groupID", "log2fold", "pvalue", "padj", "exonBaseMean"), + names(dex_df) + ), drop = FALSE] + tx_dtu <- workflow_glue_r_normalise_tsv_df(tx_dtu) + + gene_q <- per_gene_q_fn(dex_res$dxr) + gene_dtu <- data.frame( + GENEID = names(gene_q), + qval = unname(gene_q), + row.names = NULL + ) + + contrast_qc$dtu_status <- "SUCCESS" + contrast_qc$dtu_significant_transcripts <- sum(tx_dtu$padj < 0.05, na.rm = TRUE) + contrast_qc$dtu_significant_genes <- sum(gene_dtu$qval < 0.05, na.rm = TRUE) + + pdf_fn(file.path(contrast_dir, "results_dtu.pdf")) + plot_ma_fn(dex_res$dxr, cex = 0.8, alpha = 0.05) + plot_disp_fn(dex_res$dxd) + dev_off_fn() + } + + utils::write.table( + dex_df, + file = file.path(contrast_dir, "results_dexseq.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + utils::write.table( + tx_dtu[order(tx_dtu$padj), ], + file = file.path(contrast_dir, "results_dtu_transcript.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + utils::write.table( + gene_dtu[order(gene_dtu$qval), ], + file = file.path(contrast_dir, "results_dtu_gene.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + utils::write.table( + contrast_samples, + file = file.path(contrast_dir, "samples_used.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + + contrast_qc_summary <- c( + sprintf("Contrast QC Summary: %s", contrast_name), + paste(rep("=", 50), collapse = ""), + "", + "Sample Information:", + sprintf(" Target level (%s): %d samples", target_level, contrast_qc$n_target), + sprintf(" Reference level (%s): %d samples", reference_level, contrast_qc$n_reference), + sprintf(" Total samples: %d", contrast_qc$n_samples), + "", + "DGE Results:", + sprintf(" Genes tested: %d", contrast_qc$genes_tested), + sprintf(" Significant (FDR < 0.05): %d", contrast_qc$dge_significant_fdr05), + sprintf(" Significant (FDR < 0.01): %d", contrast_qc$dge_significant_fdr01), + sprintf(" Upregulated: %d", contrast_qc$dge_upregulated), + sprintf(" Downregulated: %d", contrast_qc$dge_downregulated), + "", + "DTU Results:", + sprintf(" Status: %s", contrast_qc$dtu_status), + sprintf(" Transcripts tested: %d", contrast_qc$transcripts_tested), + if (contrast_qc$dtu_status == "SUCCESS") { + c( + sprintf( + " Significant transcripts (FDR < 0.05): %d", + contrast_qc$dtu_significant_transcripts + ), + sprintf(" Genes with DTU (q < 0.05): %d", contrast_qc$dtu_significant_genes) + ) + } else { + " See DTU_ANALYSIS_FAILED.txt for details" + }, + if (!is.null(contrast_qc$dtu_power_warning)) paste0(" WARNING: ", contrast_qc$dtu_power_warning) else NULL, + "" + ) + writeLines(contrast_qc_summary, file.path(contrast_dir, "contrast_qc_summary.txt")) + de_qc_stats$contrasts[[contrast_name]] <- contrast_qc + } + + jsonlite::write_json( + de_qc_stats, + file.path(argv$out_dir, "de_qc_stats.json"), + pretty = TRUE, + auto_unbox = TRUE + ) + + overall_summary <- c( + "Differential Expression/Usage Analysis Summary", + paste(rep("=", 50), collapse = ""), + "", + sprintf("Timestamp: %s", de_qc_stats$timestamp), + sprintf("Total samples: %d", de_qc_stats$total_samples), + sprintf("Condition column: %s", de_qc_stats$condition_column), + sprintf("Reference level: %s", de_qc_stats$reference_level), + sprintf("Covariates: %s", paste(de_qc_stats$covariates, collapse = ", ")), + "", + "Sample Sizes:", + sapply(names(de_qc_stats$samples_per_group), function(grp) { + sprintf(" %s: %d samples", grp, de_qc_stats$samples_per_group[[grp]]) + }), + if (de_qc_stats$sample_size_warnings != "none") paste0(" WARNING: ", de_qc_stats$sample_size_warnings) else NULL, + "", + sprintf("Number of contrasts tested: %d", de_qc_stats$num_contrasts), + if (!is.null(de_qc_stats$multiple_testing_note)) paste0(" NOTE: ", de_qc_stats$multiple_testing_note) else NULL, + "", + "Per-Contrast Results:", + sapply(names(de_qc_stats$contrasts), function(cname) { + cqc <- de_qc_stats$contrasts[[cname]] + c( + "", + sprintf(" %s:", cname), + sprintf(" Samples: %d (%d vs %d)", cqc$n_samples, cqc$n_target, cqc$n_reference), + sprintf(" DGE significant: %d genes (FDR<0.05)", cqc$dge_significant_fdr05), + sprintf(" DTU status: %s", cqc$dtu_status), + if (cqc$dtu_status == "SUCCESS") sprintf(" DTU significant: %d genes", cqc$dtu_significant_genes) else NULL + ) + }), + "", + "For detailed per-contrast statistics, see:", + " - /contrast_qc_summary.txt", + " - /results_dge.tsv", + " - /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")) + + 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) +} diff --git a/bin/workflow_glue_r/load.R b/bin/workflow_glue_r/load.R new file mode 100644 index 0000000..3707dfc --- /dev/null +++ b/bin/workflow_glue_r/load.R @@ -0,0 +1,16 @@ +workflow_glue_r_load <- function(pkg_dir = NULL, env = globalenv()) { + if (is.null(pkg_dir)) { + stop("pkg_dir must be provided when loading workflow_glue_r.", call. = FALSE) + } + + r_dir <- file.path(pkg_dir, "R") + if (!dir.exists(r_dir)) { + stop(sprintf("R source directory not found: %s", r_dir), call. = FALSE) + } + + for (path in sort(list.files(r_dir, pattern = "\\.[Rr]$", full.names = TRUE))) { + sys.source(path, envir = env) + } + + invisible(pkg_dir) +} diff --git a/bin/workflow_glue_r/tests/testthat.R b/bin/workflow_glue_r/tests/testthat.R new file mode 100644 index 0000000..f4a710b --- /dev/null +++ b/bin/workflow_glue_r/tests/testthat.R @@ -0,0 +1,36 @@ +args <- commandArgs(trailingOnly = FALSE) +file_arg <- grep("^--file=", args, value = TRUE) +if (length(file_arg) < 1) { + stop("Unable to determine testthat.R path", call. = FALSE) +} + +test_dir <- dirname(normalizePath(sub("^--file=", "", file_arg[[1]]))) +pkg_dir <- normalizePath(file.path(test_dir, "..")) +repo_root <- normalizePath(file.path(pkg_dir, "..", "..")) + +source(file.path(pkg_dir, "load.R")) +workflow_glue_r_load(pkg_dir) + +Sys.setenv( + WORKFLOW_GLUE_R_PACKAGE_DIR = pkg_dir, + WORKFLOW_GLUE_R_REPO_ROOT = repo_root, + TEST_DATA = Sys.getenv("TEST_DATA", unset = file.path(repo_root, "test_data")) +) + +testthat_dir <- file.path(pkg_dir, "tests", "testthat") +test_files <- if (dir.exists(testthat_dir)) { + list.files(testthat_dir, pattern = "\\.[Rr]$", full.names = TRUE) +} else { + character(0) +} + +if (length(test_files) < 1) { + message("No workflow-local R test files found; skipping.") + quit(save = "no", status = 0) +} + +testthat::test_dir( + testthat_dir, + reporter = "summary", + stop_on_failure = TRUE +) diff --git a/nextflow.config b/nextflow.config index 69aaa5d..e1a6bdc 100644 --- a/nextflow.config +++ b/nextflow.config @@ -60,7 +60,7 @@ params { "--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv'", ] common_sha = "sha21d552f9910c575766e5d465fcb7b52fefda4b79" - container_sha = "shae31de64635251eca32c3f5619cb144a27f8233c6" + container_sha = "shaff0012055c9e1e71caf5b7857d717d17a3465a77" pychopper_sha = "shaaaf20a5a0e76f9e18bad21af639a6b69e4a31a2f" sqanti_sha = "sha5bd775836492699e2537ebf846098eb117191d87" agent = null diff --git a/subworkflows/differential_expression.nf b/subworkflows/differential_expression.nf index ce62c48..3caedd4 100644 --- a/subworkflows/differential_expression.nf +++ b/subworkflows/differential_expression.nf @@ -37,7 +37,7 @@ process runDifferentialAnalysis { String covariates_arg = params.covariates ? "--covariates '${params.covariates}'" : "" String reference_arg = params.reference_level ? "--reference_level '${params.reference_level}'" : "" """ - run_de_analysis.R \ + supeRglue de_analysis \ --transcript_rds "${transcript_rds}" \ --gene_rds "${gene_rds}" \ --sample_sheet "${sample_sheet}" \ diff --git a/subworkflows/transcriptome.nf b/subworkflows/transcriptome.nf index 778f5ae..908329e 100644 --- a/subworkflows/transcriptome.nf +++ b/subworkflows/transcriptome.nf @@ -179,7 +179,7 @@ process runJointBambu { String sample_sheet_arg = sample_sheet.name == OPTIONAL_FILE.name ? "" : "--sample_sheet ${sample_sheet}" String ndr_arg = params.ndr != null ? "--ndr ${params.ndr}" : "" """ - run_bambu.R \ + supeRglue bambu \ --bam_dir bams \ ${sample_sheet_arg} \ --annotation "${annotation}" \ @@ -210,7 +210,7 @@ process runPerSampleBambu { tuple val(meta), path("${meta.alias}/transcript_metadata.tsv"), emit: transcript_metadata script: """ - run_bambu.R \ + supeRglue bambu \ --bam_path "${bam}" \ --sample_alias "${meta.alias}" \ --annotation "${annotation}" \