diff --git a/bin/workflow_glue_r/R/bambu.R b/bin/workflow_glue_r/R/bambu.R index 8c22472..0c6198d 100644 --- a/bin/workflow_glue_r/R/bambu.R +++ b/bin/workflow_glue_r/R/bambu.R @@ -565,6 +565,7 @@ bambu_run_collate_mode <- function(args, gene_expression_fn, write_gtf_fn) { out_dir = args$out_dir, transcriptome_mode = args$transcriptome_mode, ndr = args$ndr, + annotation = args$annotation, gene_expression_fn = gene_expression_fn, write_gtf_fn = write_gtf_fn ) @@ -1038,6 +1039,7 @@ bambu_collate_chunk_outputs <- function( out_dir, transcriptome_mode = "discover", ndr = NULL, + annotation = NULL, gene_expression_fn = bambu::transcriptToGeneExpression, write_gtf_fn = bambu::writeToGTF ) { @@ -1077,7 +1079,8 @@ bambu_collate_chunk_outputs <- function( args <- list( out_dir = out_dir, transcriptome_mode = transcriptome_mode, - ndr = ndr + ndr = ndr, + annotation = annotation ) bambu_write_outputs( se, @@ -1111,6 +1114,62 @@ bambu_write_matrix_tsv <- function(se_obj, assay_name, id_col, meta_df, output_p invisible(gc(verbose = FALSE)) } +#' Enrich SummarizedExperiment rowData with annotation-derived names. +#' +#' Reads `gene_name` and `transcript_name` from a GFF annotation file +#' and adds them as columns in the `rowData` of the transcript-level and +#' gene-level SummarizedExperiment objects. Matching is performed on +#' `GENEID`or TXNAME +#' If `annotation` is missing/NULL, or if the GTF does not contain the +#' relevant attributes, the objects are returned unchanged. +#' +#' @param se A `SummarizedExperiment` of transcript-level counts +#' Must have GENEID in rowData for gene name enrichment. +#' @param gene_se A SummarizedExperiment of gene-level counts (bambu +#' Must have GENEID in rowData or use rownames as gene IDs. +#' @param annotation Path to a GTF/GFF annotation file, or NULL. When +#' NULL or the file is absent the function is a no-op. +#' +#' @return A named list with elements `se`` and `gene_se`, each being +#' the (possibly enriched) input object. +bambu_add_annotation_names <- function(se, gene_se, annotation) { + feature_maps <- workflow_glue_r_annotation_name_maps(annotation) + gene_name_map <- feature_maps$gene + transcript_name_map <- feature_maps$transcript + + if (nrow(gene_name_map) < 1 && nrow(transcript_name_map) < 1) { + return(list(se = se, gene_se = gene_se)) + } + + tx_row_data <- SummarizedExperiment::rowData(se) + update_tx_row_data <- FALSE + if (nrow(gene_name_map) > 0 && "GENEID" %in% names(tx_row_data)) { + tx_row_data$gene_name <- gene_name_map$gene_name[ + match(as.character(tx_row_data$GENEID), gene_name_map$GENEID) + ] + update_tx_row_data <- TRUE + } + if (nrow(transcript_name_map) > 0 && "TXNAME" %in% names(tx_row_data)) { + tx_row_data$transcript_name <- transcript_name_map$transcript_name[ + match(as.character(tx_row_data$TXNAME), transcript_name_map$TXNAME) + ] + update_tx_row_data <- TRUE + } + if (update_tx_row_data) { + SummarizedExperiment::rowData(se) <- S4Vectors::DataFrame(tx_row_data) + } + + gene_row_data <- SummarizedExperiment::rowData(gene_se) + if (nrow(gene_name_map) > 0 && "GENEID" %in% names(gene_row_data)) { + gene_row_data$gene_name <- gene_name_map$gene_name[ + match(as.character(gene_row_data$GENEID), gene_name_map$GENEID) + ] + SummarizedExperiment::rowData(gene_se) <- S4Vectors::DataFrame(gene_row_data) + } + + list(se = se, gene_se = gene_se) +} + bambu_write_outputs <- function( se, gene_se, @@ -1120,6 +1179,10 @@ bambu_write_outputs <- function( write_gtf_fn = bambu::writeToGTF, write_rds = TRUE ) { + outputs <- bambu_add_annotation_names(se, gene_se, args$annotation) + se <- outputs$se + gene_se <- outputs$gene_se + row_ranges <- SummarizedExperiment::rowRanges(se) if (length(row_ranges) > 0) { write_gtf_fn( diff --git a/bin/workflow_glue_r/R/common.R b/bin/workflow_glue_r/R/common.R index c78524e..77b334b 100644 --- a/bin/workflow_glue_r/R/common.R +++ b/bin/workflow_glue_r/R/common.R @@ -53,3 +53,79 @@ workflow_glue_r_empty_tsv <- function(columns) { names(out) <- columns out } + +#' Extract a deduplicated id-to-name mapping from a GFF annotation file. +#' +#' Uses annotation metadata object from rtracklayer to generate a two-column +#' `data.frame` mapping feature IDs to display names. +#' When a feature ID maps to multiple names only the first observed name is kept. +#' +#' @param annotation_path data.frame with annotation metadata +#' @param id_column Name of the GFF attribute to use as the identifier +#' (e.g. `gene_id` or `transcript_id`). +#' @param name_column Name of the GFF attribute to use as the display name +#' (e.g. `gene_name` or `transcript_name`). +#' @param output_id_column Column name for the identifier in the returned +#' `data.frame` (e.g. `GENEID` or `TXNAME`). +#' +#' @return A `data.frame` with columns `output_id_column` and `name_column`. +#' Returns an empty `data.frame` with those columns if the required +#' attributes are absent or all values are missing. +workflow_glue_r_annotation_name_map_from_meta <- function( + annotation_meta, + id_column, + name_column, + output_id_column +) { + output_columns <- c(output_id_column, name_column) + + if (!all(c(id_column, name_column) %in% names(annotation_meta))) { + return(workflow_glue_r_empty_tsv(output_columns)) + } + + # Keep annotation row rows only where both the id and name are present and non-empty. + keep <- + !is.na(annotation_meta[[id_column]]) & + nzchar(annotation_meta[[id_column]]) & + !is.na(annotation_meta[[name_column]]) & + nzchar(annotation_meta[[name_column]]) + if (!any(keep)) { + return(workflow_glue_r_empty_tsv(output_columns)) + } + + annotation_meta <- annotation_meta[keep, , drop = FALSE] + annotation_meta <- annotation_meta[!duplicated(annotation_meta[[id_column]]), ] + out <- data.frame( + annotation_meta[[id_column]], + annotation_meta[[name_column]], + row.names = NULL, + stringsAsFactors = FALSE + ) + names(out) <- output_columns + out +} + +#' Extract `GENEID->gene_name` and `TXNAME->transcript_name` mappings. +#' +#' @param annotation_path Path to a GTF or GFF file. +#' +#' @return A list with `gene` and `transcript` data.frames. +workflow_glue_r_annotation_name_maps <- function(annotation_path) { + if (bambu_missing(annotation_path)) { + return(list( + gene = workflow_glue_r_empty_tsv(c("GENEID", "gene_name")), + transcript = workflow_glue_r_empty_tsv(c("TXNAME", "transcript_name")) + )) + } + annotation <- rtracklayer::import(annotation_path) + annotation_meta <- S4Vectors::mcols(annotation) + + list( + gene = workflow_glue_r_annotation_name_map_from_meta( + annotation_meta, "gene_id", "gene_name", "GENEID" + ), + transcript = workflow_glue_r_annotation_name_map_from_meta( + annotation_meta, "transcript_id", "transcript_name", "TXNAME" + ) + ) +} \ No newline at end of file diff --git a/bin/workflow_glue_r/R/de_analysis.R b/bin/workflow_glue_r/R/de_analysis.R index 7b10b22..cfca013 100644 --- a/bin/workflow_glue_r/R/de_analysis.R +++ b/bin/workflow_glue_r/R/de_analysis.R @@ -672,11 +672,13 @@ de_run_dexseq_result <- function( run_inner(covariates) } -de_dge_columns <- c("GENEID", "baseMean", "log2FoldChange", "lfcSE", "stat", "pvalue", "padj") +de_dge_columns <- c("GENEID", "gene_name", "baseMean", "log2FoldChange", "lfcSE", "stat", "pvalue", "padj") de_dtu_transcript_columns <- c( "featureID", "groupID", + "gene_name", + "transcript_name", "log2FoldChange", "pvalue", "padj", @@ -1143,6 +1145,15 @@ main_run_de_analysis <- function(args) { } else { dex_df <- postprocess_run$dex_df tx_dtu <- postprocess_run$tx_dtu + if ("featureID" %in% names(tx_dtu)) { + tx_match <- match(as.character(tx_dtu$featureID), as.character(tx_meta$TXNAME)) + if ("gene_name" %in% names(tx_meta)) { + tx_dtu$gene_name <- tx_meta$gene_name[tx_match] + } + if ("transcript_name" %in% names(tx_meta)) { + tx_dtu$transcript_name <- tx_meta$transcript_name[tx_match] + } + } gene_dtu <- postprocess_run$gene_dtu contrast_qc$dtu_status <- "SUCCESS" diff --git a/bin/workflow_glue_r/tests/testthat/test_bambu.R b/bin/workflow_glue_r/tests/testthat/test_bambu.R index 47edf4c..0f654d0 100644 --- a/bin/workflow_glue_r/tests/testthat/test_bambu.R +++ b/bin/workflow_glue_r/tests/testthat/test_bambu.R @@ -768,10 +768,59 @@ testthat::test_that("list columns flattened for TSV output", { testthat::expect_equal(normalised$list_col, c("x;y", "z")) }) +testthat::test_that("annotation name maps are extracted from GTF", { + fixture_dir <- tempfile("gene-name-map-") + dir.create(fixture_dir) + + gtf <- file.path(fixture_dir, "annotation.gtf") + writeLines( + c( + paste( + "chr1", "sim", "transcript", "1", "100", ".", "+", ".", + 'gene_id "gene1"; transcript_id "tx1"; gene_name "GENEA"; transcript_name "TXA";', + sep = "\t" + ), + paste( + "chr1", "sim", "exon", "1", "100", ".", "+", ".", + 'gene_id "gene1"; transcript_id "tx1"; gene_name "GENEA"; transcript_name "TXA";', + sep = "\t" + ), + paste( + "chr1", "sim", "transcript", "201", "300", ".", "+", ".", + 'gene_id "gene2"; transcript_id "tx2"; gene_name "GENEB"; transcript_name "TXB";', + sep = "\t" + ) + ), + gtf + ) + + maps <- workflow_glue_r_annotation_name_maps(gtf) + gene_name_map <- maps$gene + transcript_name_map <- maps$transcript + + testthat::expect_equal(names(gene_name_map), c("GENEID", "gene_name")) + testthat::expect_equal(nrow(gene_name_map), 2) + testthat::expect_equal(gene_name_map$gene_name[match("gene1", gene_name_map$GENEID)], "GENEA") + testthat::expect_equal(gene_name_map$gene_name[match("gene2", gene_name_map$GENEID)], "GENEB") + + testthat::expect_equal(names(transcript_name_map), c("TXNAME", "transcript_name")) + testthat::expect_equal(nrow(transcript_name_map), 2) + testthat::expect_equal( + transcript_name_map$transcript_name[match("tx1", transcript_name_map$TXNAME)], + "TXA" + ) + testthat::expect_equal( + transcript_name_map$transcript_name[match("tx2", transcript_name_map$TXNAME)], + "TXB" + ) +}) + # 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) + fixture_dir <- tempfile("bambu-write-fixture-") + dir.create(fixture_dir) sample_names <- c("sampleA", "sampleB") base_se <- make_test_tx_se(sample_names = sample_names) @@ -782,7 +831,29 @@ testthat::test_that("bambu outputs written correctly", { ) gene_se <- make_test_gene_se(sample_names = sample_names) sample_df <- data.frame(alias = sample_names, stringsAsFactors = FALSE) - args <- list(out_dir = out_dir, transcriptome_mode = "discover", ndr = 0.15) + annotation <- file.path(fixture_dir, "gene_names.gtf") + writeLines( + c( + paste( + "chr1", "test", "transcript", "1", "50", ".", "+", ".", + 'gene_id "gene1"; transcript_id "tx1"; gene_name "GENEA"; transcript_name "TXA";', + sep = "\t" + ), + paste( + "chr1", "test", "transcript", "201", "250", ".", "+", ".", + 'gene_id "gene2"; transcript_id "tx3"; gene_name "GENEB"; transcript_name "TXC";', + sep = "\t" + ) + ), + annotation + ) + + args <- list( + out_dir = out_dir, + transcriptome_mode = "discover", + ndr = 0.15, + annotation = annotation + ) qc_stats <- list( samples = 2, total_transcripts_before_filter = 4, @@ -828,7 +899,37 @@ testthat::test_that("bambu outputs written correctly", { testthat::expect_true("eqClassById" %in% names(tx_meta)) testthat::expect_equal(tx_meta$eqClassById[[1]], "1;2") + testthat::expect_true("gene_name" %in% names(tx_meta)) + testthat::expect_true("transcript_name" %in% names(tx_meta)) + testthat::expect_equal( + unique(stats::na.omit(tx_meta$gene_name[tx_meta$GENEID == "gene1"])), + "GENEA" + ) + testthat::expect_equal( + tx_meta$transcript_name[match("tx1", tx_meta$TXNAME)], + "TXA" + ) testthat::expect_true(all(c("TXNAME", "sampleA", "sampleB") %in% names(tx_counts))) + + tx_rds <- readRDS(file.path(out_dir, "bambu_transcripts.rds")) + gene_rds <- readRDS(file.path(out_dir, "bambu_genes.rds")) + tx_rds_meta <- as.data.frame(SummarizedExperiment::rowData(tx_rds)) + gene_rds_meta <- as.data.frame(SummarizedExperiment::rowData(gene_rds)) + testthat::expect_true("gene_name" %in% names(tx_rds_meta)) + testthat::expect_true("transcript_name" %in% names(tx_rds_meta)) + testthat::expect_true("gene_name" %in% names(gene_rds_meta)) + testthat::expect_equal( + unique(stats::na.omit(tx_rds_meta$gene_name[tx_rds_meta$GENEID == "gene2"])), + "GENEB" + ) + testthat::expect_equal( + tx_rds_meta$transcript_name[match("tx3", tx_rds_meta$TXNAME)], + "TXC" + ) + testthat::expect_equal( + gene_rds_meta$gene_name[match("gene1", gene_rds_meta$GENEID)], + "GENEA" + ) }) testthat::test_that("collate combines multiple chunk quantification outputs", { diff --git a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R b/bin/workflow_glue_r/tests/testthat/test_de_analysis.R index b404eec..e78e625 100644 --- a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R +++ b/bin/workflow_glue_r/tests/testthat/test_de_analysis.R @@ -813,6 +813,19 @@ testthat::test_that("CLI integration produces expected outputs", { fixture_dir <- tempfile("de-cli-") dir.create(fixture_dir) bundle <- write_de_fixture_bundle(fixture_dir, levels = c("control", "treated")) + tx_se <- readRDS(bundle$transcript_rds) + tx_row_data <- as.data.frame(SummarizedExperiment::rowData(tx_se)) + tx_row_data$gene_name <- c("GENEA", "GENEA", "GENEB", "GENEB") + tx_row_data$transcript_name <- c("TXA", "TXB", "TXC", "TXD") + SummarizedExperiment::rowData(tx_se) <- S4Vectors::DataFrame(tx_row_data) + saveRDS(tx_se, bundle$transcript_rds) + + gene_se <- readRDS(bundle$gene_rds) + gene_row_data <- as.data.frame(SummarizedExperiment::rowData(gene_se)) + gene_row_data$gene_name <- c("GENEA", "GENEB") + SummarizedExperiment::rowData(gene_se) <- S4Vectors::DataFrame(gene_row_data) + saveRDS(gene_se, bundle$gene_rds) + out_dir <- file.path(fixture_dir, "out") result <- run_rscript( @@ -850,8 +863,11 @@ testthat::test_that("CLI integration produces expected outputs", { ) testthat::expect_gt(nrow(dge), 0) - testthat::expect_true(all(c("GENEID", "log2FoldChange", "padj") %in% names(dge))) - testthat::expect_true(all(c("featureID", "groupID", "padj") %in% names(dtu_tx))) + testthat::expect_true(all(c("GENEID", "gene_name", "log2FoldChange", "padj") %in% names(dge))) + testthat::expect_true(all(c("featureID", "groupID", "gene_name", "transcript_name", "padj") %in% names(dtu_tx))) + testthat::expect_true(all(stats::na.omit(dge$gene_name) %in% c("GENEA", "GENEB"))) + testthat::expect_true(all(stats::na.omit(dtu_tx$gene_name) %in% c("GENEA", "GENEB"))) + testthat::expect_true(all(stats::na.omit(dtu_tx$transcript_name) %in% c("TXA", "TXB", "TXC", "TXD"))) testthat::expect_true(all(c("featureID", "groupID", "padj") %in% names(dexseq))) testthat::expect_true("analysis_fallbacks" %in% names(de_qc)) testthat::expect_true("contrasts" %in% names(de_qc)) diff --git a/modules/local/bambu_chunked.nf b/modules/local/bambu_chunked.nf index 92bcd1e..aa49e47 100644 --- a/modules/local/bambu_chunked.nf +++ b/modules/local/bambu_chunked.nf @@ -88,6 +88,7 @@ process collateBambuQuant { memory "16 GB" input: tuple val(meta), path(chunk_dirs, stageAs: "chunks/*") + path annotation, stageAs: "annotation/*" output: tuple val(meta), path("${meta.alias}"), emit: dir tuple val(meta), path("${meta.alias}/transcripts.gtf"), emit: gtf @@ -103,6 +104,7 @@ process collateBambuQuant { """ supeRglue bambu collate \ ${chunk_dirs_arg} \ + --annotation "${annotation}" \ --transcriptome_mode "${params.transcriptome_mode}" \ ${ndr_arg} \ --out_dir "${meta.alias}" diff --git a/subworkflows/transcriptome.nf b/subworkflows/transcriptome.nf index d139918..c41d45b 100644 --- a/subworkflows/transcriptome.nf +++ b/subworkflows/transcriptome.nf @@ -252,7 +252,8 @@ workflow transcriptome_analysis { .groupTuple() .map { alias, metas, chunk_dirs -> tuple(metas[0], chunk_dirs) - } + }, + analysis_annotation ) joint_bambu_empty = runJointBambuEmpty( bambu_empty_inputs(joint_quant_inputs_all) @@ -294,7 +295,8 @@ workflow transcriptome_analysis { .groupTuple() .map { alias, metas, chunk_dirs -> tuple(metas[0], chunk_dirs) - } + }, + analysis_annotation ) sample_bambu_empty = runPerSampleBambuEmpty( bambu_empty_inputs(sample_quant_inputs_all)