Merge branch 'chunk-quant' into 'dev'

Scatter-gather bambu quant [CW-7201]

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

bin/workflow_glue/report.py

@@ -1,6 +1,7 @@
 """Create workflow report for wf-transcriptomes."""
 
 import json
+import math
 from pathlib import Path
 
 from dominate.tags import div, h3, p, pre, strong
@@ -25,6 +26,35 @@ def _read_table(path, **kwargs):
     return pd.read_csv(path, sep="\t", **kwargs)
 
 
+def _coerce_float(value):
+    """Return a finite float when possible, otherwise None."""
+    if value in (None, "", "N/A", "NA", "nan", "NaN"):
+        return None
+    try:
+        numeric = float(value)
+    except (TypeError, ValueError):
+        return None
+    if not math.isfinite(numeric):
+        return None
+    return numeric
+
+
+def _format_count_value(value):
+    """Format count-like values for the report, tolerating NA-like strings."""
+    numeric = _coerce_float(value)
+    if numeric is None:
+        return "N/A"
+    return format(round(numeric), ",")
+
+
+def _format_ratio_value(value):
+    """Format ratio values for the report, tolerating NA-like strings."""
+    numeric = _coerce_float(value)
+    if numeric is None:
+        return "N/A"
+    return f"{numeric:.2f}x"
+
+
 def _cohort_summary(cohort_dir):
     tx_meta = _read_table(Path(cohort_dir) / "transcript_metadata.tsv")
     if tx_meta is None:
@@ -262,12 +292,23 @@ def main(args):
                 stats = stats[0]
                 flagstats = flagstats[0]
                 sample_names = sample_names[0] if sample_names else None
-            fastcat.SeqSummary(
-                stats,
-                flagstat=flagstats,
-                sample_names=sample_names,
-                alignment_stats=True,
-            )
+            try:
+                fastcat.SeqSummary(
+                    stats,
+                    flagstat=flagstats,
+                    sample_names=sample_names,
+                    alignment_stats=True,
+                )
+            except Exception as exc:  # pragma: no cover - defensive
+                logger.warning("Skipping read summary plot: %s", exc)
+                _create_warning_banner(
+                    (
+                        "Read summary plots could not be rendered for this run. "
+                        "This can happen for degenerate or extremely "
+                        "small input statistics."
+                    ),
+                    level="info",
+                )
 
     with report.add_section("Sample metadata", "Samples"):
         tabs = Tabs()
@@ -389,33 +430,30 @@ def main(args):
                         (
                             "Median library size",
                             "{} reads".format(
-                                format(
-                                    bambu_qc.get("median_library_size", 0),
-                                    ",",
+                                _format_count_value(
+                                    bambu_qc.get("median_library_size", 0)
                                 )
                             ),
                         ),
                         (
                             "Min library size",
                             "{} reads".format(
-                                format(
-                                    bambu_qc.get("min_library_size", 0),
-                                    ",",
+                                _format_count_value(
+                                    bambu_qc.get("min_library_size", 0)
                                 )
                             ),
                         ),
                         (
                             "Max library size",
                             "{} reads".format(
-                                format(
-                                    bambu_qc.get("max_library_size", 0),
-                                    ",",
+                                _format_count_value(
+                                    bambu_qc.get("max_library_size", 0)
                                 )
                             ),
                         ),
                         (
                             "Library size ratio (max/min)",
-                            "{:.2f}x".format(
+                            _format_ratio_value(
                                 bambu_qc.get("library_size_ratio", 1.0)
                             ),
                         ),
@@ -463,11 +501,14 @@ def main(args):
                 with h3("Per-Sample Library Sizes"):
                     lib_size_data = []
                     for sample, size in bambu_qc["library_sizes"].items():
+                        numeric_size = _coerce_float(size)
                         lib_size_data.append(
                             {
                                 "Sample": sample,
-                                "Library Size": format(size, ","),
-                                "Reads": size,
+                                "Library Size": _format_count_value(size),
+                                "Reads": (
+                                    numeric_size if numeric_size is not None else -1
+                                ),
                             }
                         )
                     lib_df = pd.DataFrame(lib_size_data).sort_values(

bin/workflow_glue/tests/common/test_report.py

@@ -217,6 +217,132 @@ def test_report_main_accepts_optional_file_sentinels(monkeypatch, tmp_path):
     assert any("GRCh38" in table.to_string() for table in tables)
 
 
+def test_report_main_handles_degenerate_bambu_qc_and_read_summary(
+    monkeypatch,
+    tmp_path,
+):
+    """Tiny/empty stats should not crash the report rendering path."""
+    tables = []
+    banners = []
+
+    monkeypatch.setattr(report.labs, "LabsReport", _FakeReport)
+    monkeypatch.setattr(report, "Tabs", _FakeTabs)
+    monkeypatch.setattr(report, "p", lambda *args, **kwargs: None)
+    monkeypatch.setattr(report, "pre", lambda *args, **kwargs: None)
+    monkeypatch.setattr(
+        report.fastcat,
+        "SeqSummary",
+        lambda *args, **kwargs: (_ for _ in ()).throw(KeyError(1)),
+    )
+    monkeypatch.setattr(
+        report,
+        "_create_warning_banner",
+        lambda message, level="warning": banners.append((level, message)),
+    )
+    monkeypatch.setattr(
+        report.DataTable,
+        "from_pandas",
+        staticmethod(lambda table, *args, **kwargs: tables.append(table.copy())),
+    )
+
+    metadata = _write(
+        tmp_path / "metadata.json",
+        json.dumps([{"alias": "sampleA", "has_stats": True}]),
+    )
+    params = _write(tmp_path / "params.json", "{}")
+    versions = tmp_path / "versions"
+    versions.mkdir()
+    _write(versions / "versions.txt", "tool,1.0\n")
+
+    cohort = tmp_path / "cohort"
+    cohort.mkdir()
+    reference = cohort / "reference"
+    reference.mkdir()
+    _write(
+        reference / "annotation_reference_summary.json",
+        json.dumps(
+            {
+                "seqname_overlap": ["chr1"],
+                "only_in_annotation": [],
+                "only_in_reference": [],
+                "annotation": {
+                    "kept_records": 10,
+                    "excluded_unstranded_records": 0,
+                    "sanitised_attribute_records": 0,
+                },
+                "warnings": [],
+            }
+        ),
+    )
+    _write(
+        cohort / "bambu_qc_stats.json",
+        json.dumps(
+            {
+                "samples": 1,
+                "library_sizes": {"sampleA": 0},
+                "min_library_size": 0,
+                "max_library_size": 0,
+                "median_library_size": 0,
+                "library_size_ratio": "NA",
+                "total_transcripts_before_filter": 0,
+                "total_transcripts_after_filter": 0,
+                "transcripts_filtered": 0,
+                "median_transcripts_detected": 0,
+                "total_genes_after_filter": 0,
+                "transcriptome_mode": "discover",
+                "ndr_used": "automatic",
+            }
+        ),
+    )
+
+    samples = tmp_path / "samples"
+    samples.mkdir()
+    (samples / "OPTIONAL_FILE").touch()
+
+    sqanti = tmp_path / "sqanti"
+    sqanti.mkdir()
+    (sqanti / "OPTIONAL_FILE").touch()
+
+    alignment_stats = tmp_path / "alignment_stats"
+    alignment_stats.mkdir()
+
+    out_report = tmp_path / "wf-transcriptomes-report.html"
+    args = report.argparser().parse_args(
+        [
+            str(out_report),
+            "--metadata",
+            str(metadata),
+            "--alignment_stats_dir",
+            str(alignment_stats),
+            "--stats",
+            str(alignment_stats),
+            "--cohort_dir",
+            str(cohort),
+            "--samples_dir",
+            str(samples),
+            "--sqanti_dir",
+            str(sqanti),
+            "--versions",
+            str(versions),
+            "--params",
+            str(params),
+        ]
+    )
+
+    report.main(args)
+
+    assert out_report.exists()
+    assert any(
+        level == "info" and "Read summary plots could not be rendered" in message
+        for level, message in banners
+    )
+    assert any(
+        "Library size ratio (max/min)" in table.to_string()
+        and "N/A" in table.to_string()
+        for table in tables
+    )
+
+
 def test_report_main_renders_statistical_methods_and_warnings(
     monkeypatch,
     tmp_path,

modules/local/bambu_chunked.nf

@@ -0,0 +1,120 @@
+nextflow.enable.dsl = 2
+
+OPTIONAL_FILE = file("$projectDir/data/OPTIONAL_FILE")
+
+
+process bambuDiscover {
+    label "wf_transcriptomes"
+    cpus {
+        int requested = (params.threads ?: 4) as int
+        int sampleCount = aliases instanceof Collection ? aliases.size() : 1
+        sampleCount > 1 ? requested : 1
+    }
+    memory "60 GB"
+    input:
+        tuple val(meta), val(aliases), path(bams, stageAs: "bams/??.bam"), path(bais, stageAs: "bams/??.bam.bai"), path(sample_sheet)
+        path annotation, stageAs: "annotation/*"
+        path reference, stageAs: "reference/*"
+    output:
+        tuple val(meta), path("discover"), emit: dir
+    script:
+        def bam_list = bams instanceof Collection ? bams : [bams]
+        def alias_list = aliases instanceof Collection ? aliases : [aliases]
+        String bams_arg = "--bams '${bam_list.join(",")}'"
+        String aliases_arg = "--aliases '${alias_list.join(",")}'"
+        String sample_sheet_arg = sample_sheet.name == OPTIONAL_FILE.name ? "" : "--sample_sheet '${sample_sheet}'"
+        String ndr_arg = params.ndr != null ? "--ndr ${params.ndr}" : ""
+    """
+    supeRglue bambu discover \
+        ${bams_arg} \
+        ${aliases_arg} \
+        ${sample_sheet_arg} \
+        --annotation "${annotation}" \
+        --genome "${reference}" \
+        --transcriptome_mode "${params.transcriptome_mode}" \
+        --threads ${task.cpus} \
+        ${ndr_arg} \
+        --out_dir discover
+    """
+}
+
+
+process bambuQuant {
+    label "wf_transcriptomes"
+    cpus {
+        int requested = (params.threads ?: 4) as int
+        boolean isJoint = meta instanceof Map && meta.alias == 'cohort'
+        isJoint ? requested : 1
+    }
+    memory { ["8.GB", "16.GB", "48.GB"][task.attempt - 1] }
+    maxRetries 2
+    errorStrategy 'retry'
+    input:
+        tuple val(meta), val(chunk_id), val(annotation_tx_count), path(chunk_rds), path(discovered_annotation)
+        path reference, stageAs: "reference/*"
+    output:
+        tuple val(meta), val(chunk_id), path("${chunk_id}"), emit: dir
+    script:
+    """
+    supeRglue bambu quant \
+        --chunk_rds "${chunk_rds}" \
+        --discovered_annotation_rds "${discovered_annotation}" \
+        --genome "${reference}" \
+        --threads ${task.cpus} \
+        --out_dir "${chunk_id}"
+    """
+}
+
+
+process bambuEmpty {
+    label "wf_transcriptomes"
+    cpus 1
+    memory "4 GB"
+    input:
+        tuple val(meta), val(aliases)
+    output:
+        tuple val(meta), path("${meta.alias}"), emit: dir
+        tuple val(meta), path("${meta.alias}/transcripts.gtf"), emit: gtf
+        tuple val(meta), path("${meta.alias}/transcript_counts.tsv"), emit: transcript_counts
+        tuple val(meta), path("${meta.alias}/gene_counts.tsv"), emit: gene_counts
+        tuple val(meta), path("${meta.alias}/bambu_transcripts.rds"), emit: transcript_rds
+        tuple val(meta), path("${meta.alias}/bambu_genes.rds"), emit: gene_rds
+        tuple val(meta), path("${meta.alias}/transcript_metadata.tsv"), emit: transcript_metadata
+    script:
+        def alias_list = aliases instanceof Collection ? aliases : [aliases]
+        String aliases_arg = "--aliases '${alias_list.join(",")}'"
+    """
+    supeRglue bambu empty \
+        ${aliases_arg} \
+        --transcriptome_mode "${params.transcriptome_mode}" \
+        --out_dir "${meta.alias}"
+    """
+}
+
+
+process collateBambuQuant {
+    label "wf_transcriptomes"
+    cpus 1
+    memory "16 GB"
+    input:
+        tuple val(meta), path(chunk_dirs, stageAs: "chunks/*")
+    output:
+        tuple val(meta), path("${meta.alias}"), emit: dir
+        tuple val(meta), path("${meta.alias}/transcripts.gtf"), emit: gtf
+        tuple val(meta), path("${meta.alias}/transcript_counts.tsv"), emit: transcript_counts
+        tuple val(meta), path("${meta.alias}/gene_counts.tsv"), emit: gene_counts
+        tuple val(meta), path("${meta.alias}/bambu_transcripts.rds"), emit: transcript_rds
+        tuple val(meta), path("${meta.alias}/bambu_genes.rds"), emit: gene_rds
+        tuple val(meta), path("${meta.alias}/transcript_metadata.tsv"), emit: transcript_metadata
+    script:
+        def chunk_dir_list = chunk_dirs instanceof Collection ? chunk_dirs : [chunk_dirs]
+        String chunk_dirs_arg = "--chunk_dirs '${chunk_dir_list.join(",")}'"
+        String ndr_arg = params.ndr != null ? "--ndr ${params.ndr}" : ""
+    """
+    supeRglue bambu collate \
+        ${chunk_dirs_arg} \
+        --transcriptome_mode "${params.transcriptome_mode}" \
+        ${ndr_arg} \
+        --out_dir "${meta.alias}"
+    """
+}

subworkflows/transcriptome.nf

@@ -2,6 +2,84 @@ nextflow.enable.dsl = 2
 
 OPTIONAL_FILE = file("$projectDir/data/OPTIONAL_FILE")
 
+// use nextflows classic rename trick
+include {
+    // joint
+    bambuDiscover as runJointBambuDiscover
+    bambuQuant as runJointBambuQuant
+    bambuEmpty as runJointBambuEmpty
+    collateBambuQuant as collateJointBambuQuant
+    // persample
+    bambuDiscover as runPerSampleBambuDiscover
+    bambuQuant as runPerSampleBambuQuant
+    bambuEmpty as runPerSampleBambuEmpty
+    collateBambuQuant as collatePerSampleBambuQuant
+} from '../modules/local/bambu_chunked'
+
+
+def bambu_discover_to_quant_inputs(discover_channel) {
+    discover_channel
+        .map { meta, discover_dir ->
+            def sampleAliases = discover_dir.resolve("samples.csv")
+                .readLines()
+                .drop(1)
+                .findAll { it?.trim() }
+                .collect { it.split(",", 2)[0] }
+            tuple(
+                meta,
+                sampleAliases,
+                discover_dir.resolve("bambu_discovered_annotations.rds"),
+                discover_dir.resolve("chunks"),
+                discover_dir.resolve("chunk_manifest.tsv")
+            )
+        }
+        .splitCsv(header: true, sep: '\t', elem: 4)
+        .map { meta, sample_aliases, discovered_annotation, chunks_dir, row ->
+            def txCountRaw = row.annotation_tx_count?.toString()?.trim()
+            Integer annotationTxCount = (!txCountRaw || txCountRaw == 'NA') ? null : txCountRaw as Integer
+            tuple(
+                meta,
+                sample_aliases,
+                row.chunk_id,
+                annotationTxCount,
+                chunks_dir.resolve("${row.chunk_id}.rds"),
+                discovered_annotation
+            )
+        }
+}
+
+
+def bambu_filter_quant_inputs_with_warning(quant_inputs) {
+    quant_inputs.filter { meta, sample_aliases, chunk_id, annotation_tx_count, chunk_rds, discovered_annotation ->
+        boolean keep = annotation_tx_count == null || annotation_tx_count > 0
+        if (!keep) {
+            log.warn("Dropping bambu quant chunk '${chunk_id}' for '${meta.alias}' because annotation_tx_count=0")
+        }
+        keep
+    }
+}
+
+
+def bambu_empty_inputs(quant_inputs) {
+    quant_inputs
+        .map { meta, sample_aliases, chunk_id, annotation_tx_count, chunk_rds, discovered_annotation ->
+            tuple(meta.alias, meta, sample_aliases, annotation_tx_count)
+        }
+        .groupTuple()
+        .filter { alias, metas, sample_aliases_sets, annotation_tx_counts ->
+            !annotation_tx_counts.any { it == null || it > 0 }
+        }
+        .map { alias, metas, sample_aliases_sets, annotation_tx_counts ->
+            tuple(metas[0], sample_aliases_sets[0])
+        }
+}
+
+
+def bambu_quant_process_inputs(quant_inputs) {
+    quant_inputs.map { meta, sample_aliases, chunk_id, annotation_tx_count, chunk_rds, discovered_annotation ->
+        tuple(meta, chunk_id, annotation_tx_count, chunk_rds, discovered_annotation)
+    }
+}
 
 process prepareAnnotationReference {
     label "wf_transcriptomes"
@@ -27,79 +105,6 @@ process prepareAnnotationReference {
 }
 
 
-process runJointBambu {
-    label "wf_transcriptomes"
-    cpus { params.threads ?: 4 }
-    memory "32 GB"
-    input:
-        tuple val(aliases), path(bams, stageAs: "bams/??.bam"), path(bais, stageAs: "bams/??.bam.bai")
-        path sample_sheet
-        path annotation, stageAs: "annotation/*"
-        path reference, stageAs: "reference/*"
-    output:
-        path "cohort", emit: dir
-        path "cohort/transcripts.gtf", emit: gtf
-        path "cohort/transcript_counts.tsv", emit: transcript_counts
-        path "cohort/gene_counts.tsv", emit: gene_counts
-        path "cohort/bambu_transcripts.rds", emit: transcript_rds
-        path "cohort/bambu_genes.rds", emit: gene_rds
-        path "cohort/transcript_metadata.tsv", emit: transcript_metadata
-    script:
-        def bam_list = bams instanceof Collection ? bams : [bams]  // todo dont run joint on single sample anyway
-        def alias_list = aliases instanceof Collection ? aliases : [aliases]
-        String bams_arg = "--bams '${bam_list.join(",")}'"
-        String aliases_arg = "--aliases '${alias_list.join(",")}'"
-        String sample_sheet_arg = sample_sheet.name == OPTIONAL_FILE.name ? "" : "--sample_sheet ${sample_sheet}"
-        String ndr_arg = params.ndr != null ? "--ndr ${params.ndr}" : ""
-    """
-    supeRglue bambu \
-        ${bams_arg} \
-        ${aliases_arg} \
-        ${sample_sheet_arg} \
-        --annotation "${annotation}" \
-        --genome "${reference}" \
-        --transcriptome_mode "${params.transcriptome_mode}" \
-        --threads ${task.cpus} \
-        ${ndr_arg} \
-        --out_dir cohort \
-    """
-}
-
-
-process runPerSampleBambu {
-    label "wf_transcriptomes"
-    cpus { params.threads ?: 4 }
-    memory "24 GB"
-    input:
-        tuple val(meta), path(bam), path(bai), path(stats)
-        path annotation, stageAs: "annotation/*"
-        path reference, stageAs: "reference/*"
-    output:
-        tuple val(meta), path("${meta.alias}"), emit: dir
-        tuple val(meta), path("${meta.alias}/transcripts.gtf"), emit: gtf
-        tuple val(meta), path("${meta.alias}/transcript_counts.tsv"), emit: transcript_counts
-        tuple val(meta), path("${meta.alias}/gene_counts.tsv"), emit: gene_counts
-        tuple val(meta), path("${meta.alias}/bambu_transcripts.rds"), emit: transcript_rds
-        tuple val(meta), path("${meta.alias}/bambu_genes.rds"), emit: gene_rds
-        tuple val(meta), path("${meta.alias}/transcript_metadata.tsv"), emit: transcript_metadata
-    script:
-        String bams_arg = "--bams '${bam.toString()}'"
-        String aliases_arg = "--aliases '${meta.alias}'"
-        String ndr_arg = params.ndr != null ? "--ndr ${params.ndr}" : ""
-    """
-    supeRglue bambu \
-        ${bams_arg} \
-        ${aliases_arg} \
-        --annotation "${annotation}" \
-        --genome "${reference}" \
-        --transcriptome_mode "${params.transcriptome_mode}" \
-        --threads ${task.cpus} \
-        ${ndr_arg} \
-        --out_dir "${meta.alias}"
-    """
-}
-
-
 process buildCohortTranscriptomeFasta {
     label "wf_transcriptomes"
     cpus 1
@@ -214,33 +219,112 @@ workflow transcriptome_analysis {
         analysis_annotation = prepared_reference_annotation.annotation.first()
         analysis_reference = prepared_reference_annotation.reference.first()
 
-        joint_bambu = runJointBambu(
+        joint_meta = [alias: "cohort"]
+
+        joint_discover = runJointBambuDiscover(
             alignments
-            | collect(flat: false)
-            | map { rows ->
-                // transform [meta, bam, bai] to [[alias1...aliasN], [bam1...bamN], [bai1...baiN]]
-                tuple(
-                    rows.collect { it[0].alias },
-                    rows.collect { it[1] },
-                    rows.collect { it[2] }
-                )
-            },
-            sample_sheet,
+                .collect(flat: false)
+                .map { rows ->
+                    // transform [meta, bam, bai] rows to 
+                    // [meta, [alias1...aliasN], [bam1...bamN], [bai1...baiN], sample_sheet]
+                    tuple(
+                        joint_meta,
+                        rows.collect { it[0].alias },
+                        rows.collect { it[1] },
+                        rows.collect { it[2] },
+                        sample_sheet
+                    )
+                },
             analysis_annotation,
             analysis_reference
         )
+        joint_quant_inputs_all = bambu_discover_to_quant_inputs(joint_discover.dir)
+        joint_quant = runJointBambuQuant(
+            bambu_quant_process_inputs(
+                bambu_filter_quant_inputs_with_warning(joint_quant_inputs_all)
+            ),
+            analysis_reference
+        )
+        joint_bambu_real = collateJointBambuQuant(
+            joint_quant.dir
+                .map { meta, chunk_id, chunk_dir ->
+                    tuple(meta.alias, meta, chunk_dir)
+                }
+                .groupTuple()
+                .map { alias, metas, chunk_dirs ->
+                    tuple(metas[0], chunk_dirs)
+                }
+        )
+        joint_bambu_empty = runJointBambuEmpty(
+            bambu_empty_inputs(joint_quant_inputs_all)
+        )
 
-        sample_bambu = runPerSampleBambu(alignments, analysis_annotation, analysis_reference)
+        joint_bambu_dir = joint_bambu_real.dir.mix(joint_bambu_empty.dir)
+        joint_bambu_gtf = joint_bambu_real.gtf.mix(joint_bambu_empty.gtf)
+        joint_bambu_transcript_counts = joint_bambu_real.transcript_counts.mix(joint_bambu_empty.transcript_counts)
+        joint_bambu_gene_counts = joint_bambu_real.gene_counts.mix(joint_bambu_empty.gene_counts)
+        joint_bambu_transcript_rds = joint_bambu_real.transcript_rds.mix(joint_bambu_empty.transcript_rds)
+        joint_bambu_gene_rds = joint_bambu_real.gene_rds.mix(joint_bambu_empty.gene_rds)
+        joint_bambu_metadata = joint_bambu_real.transcript_metadata.mix(joint_bambu_empty.transcript_metadata)
 
-        joint_fasta = buildCohortTranscriptomeFasta(joint_bambu.gtf, analysis_reference)
-        sample_fastas = buildSampleTranscriptomeFasta(sample_bambu.gtf, analysis_reference)
+        sample_discover = runPerSampleBambuDiscover(
+            alignments.map { meta, bam, bai, stats ->
+                tuple(
+                    meta,
+                    [meta.alias],
+                    bam,
+                    bai,
+                    sample_sheet
+                )
+            },
+            analysis_annotation,
+            analysis_reference
+        )
+        sample_quant_inputs_all = bambu_discover_to_quant_inputs(sample_discover.dir)
+        sample_quant = runPerSampleBambuQuant(
+            bambu_quant_process_inputs(
+                bambu_filter_quant_inputs_with_warning(sample_quant_inputs_all)
+            ),
+            analysis_reference
+        )
+        sample_bambu_real = collatePerSampleBambuQuant(
+            sample_quant.dir
+                .map { meta, chunk_id, chunk_dir ->
+                    tuple(meta.alias, meta, chunk_dir)
+                }
+                .groupTuple()
+                .map { alias, metas, chunk_dirs ->
+                    tuple(metas[0], chunk_dirs)
+                }
+        )
+        sample_bambu_empty = runPerSampleBambuEmpty(
+            bambu_empty_inputs(sample_quant_inputs_all)
+        )
+
+        sample_bambu_dirs = sample_bambu_real.dir.mix(sample_bambu_empty.dir)
+        sample_bambu_gtf = sample_bambu_real.gtf.mix(sample_bambu_empty.gtf)
+        sample_bambu_transcript_counts = sample_bambu_real.transcript_counts.mix(sample_bambu_empty.transcript_counts)
+        sample_bambu_gene_counts = sample_bambu_real.gene_counts.mix(sample_bambu_empty.gene_counts)
+        sample_bambu_transcript_rds = sample_bambu_real.transcript_rds.mix(sample_bambu_empty.transcript_rds)
+        sample_bambu_gene_rds = sample_bambu_real.gene_rds.mix(sample_bambu_empty.gene_rds)
+        sample_bambu_metadata = sample_bambu_real.transcript_metadata.mix(sample_bambu_empty.transcript_metadata)
+
+        joint_fasta = buildCohortTranscriptomeFasta(
+            joint_bambu_real.gtf.map { meta, gtf -> gtf },
+            analysis_reference
+        )
+        sample_fastas = buildSampleTranscriptomeFasta(sample_bambu_real.gtf, analysis_reference)
 
         if (params.skip_sqanti) {
             joint_sqanti_dir = Channel.empty()
             sample_sqanti_dirs = Channel.empty()
         } else {
-            joint_sqanti = runJointSqanti(joint_bambu.gtf, analysis_annotation, analysis_reference)
-            sample_sqanti = runPerSampleSqanti(sample_bambu.gtf, analysis_annotation, analysis_reference)
+            joint_sqanti = runJointSqanti(
+                joint_bambu_real.gtf.map { meta, gtf -> gtf },
+                analysis_annotation,
+                analysis_reference
+            )
+            sample_sqanti = runPerSampleSqanti(sample_bambu_real.gtf, analysis_annotation, analysis_reference)
             joint_sqanti_dir = joint_sqanti.dir
             sample_sqanti_dirs = sample_sqanti.dir
         }
@@ -249,22 +333,22 @@ workflow transcriptome_analysis {
         annotation = analysis_annotation
         annotation_reference_summary = prepared_reference_annotation.summary
         unstranded_annotation = prepared_reference_annotation.unstranded
-        joint_dir = joint_bambu.dir
-        joint_gtf = joint_bambu.gtf
+        joint_dir = joint_bambu_dir.map { meta, dir -> dir }
+        joint_gtf = joint_bambu_gtf.map { meta, gtf -> gtf }
         joint_fasta = joint_fasta.fasta
-        joint_transcript_counts = joint_bambu.transcript_counts
-        joint_gene_counts = joint_bambu.gene_counts
-        joint_transcript_rds = joint_bambu.transcript_rds
-        joint_gene_rds = joint_bambu.gene_rds
-        joint_metadata = joint_bambu.transcript_metadata
-        sample_dirs = sample_bambu.dir
-        sample_gtf = sample_bambu.gtf
+        joint_transcript_counts = joint_bambu_transcript_counts.map { meta, counts -> counts }
+        joint_gene_counts = joint_bambu_gene_counts.map { meta, counts -> counts }
+        joint_transcript_rds = joint_bambu_transcript_rds.map { meta, rds -> rds }
+        joint_gene_rds = joint_bambu_gene_rds.map { meta, rds -> rds }
+        joint_metadata = joint_bambu_metadata.map { meta, metadata -> metadata }
+        sample_dirs = sample_bambu_dirs
+        sample_gtf = sample_bambu_gtf
         sample_fastas = sample_fastas.fasta
-        sample_transcript_counts = sample_bambu.transcript_counts
-        sample_gene_counts = sample_bambu.gene_counts
-        sample_transcript_rds = sample_bambu.transcript_rds
-        sample_gene_rds = sample_bambu.gene_rds
-        sample_metadata = sample_bambu.transcript_metadata
+        sample_transcript_counts = sample_bambu_transcript_counts
+        sample_gene_counts = sample_bambu_gene_counts
+        sample_transcript_rds = sample_bambu_transcript_rds
+        sample_gene_rds = sample_bambu_gene_rds
+        sample_metadata = sample_bambu_metadata
         joint_sqanti_dir = joint_sqanti_dir
         sample_sqanti_dirs = sample_sqanti_dirs
 }