Merge branch 'gene_names' into 'dev'

Add gene name to TSV outputs [CW-7160]

See merge request epi2melabs/workflows/wf-transcriptomes!274

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(

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"
+        )
+    )
+}

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"

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", {

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("GENEID", "gene_name", "log2FoldChange", "padj") %in% names(dge)))
-    testthat::expect_true(all(c("featureID", "groupID", "padj") %in% names(dtu_tx)))
+    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))

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

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)