Use new alignment ingress [CW-7186][CW-7202][CW-7229]

2026-05-14 09:12:38 +00:00 · 2026-05-14 09:12:38 +00:00 · 06d77b2543
commit 06d77b2543
parent a99e5cda04
23 changed files with 942 additions and 658 deletions
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@ -15,6 +15,7 @@ variables:
              PYTEST_CONTAINER_CONFIG_KEY: "common_sha"
              RTEST_CONTAINER_NAME: "wf-transcriptomes-core"
              RTEST_CONTAINER_CONFIG_KEY: "container_sha"
              WF_TEMPLATE_ENFORCEMENT_BRANCH: "CW-6552"
 macos-run:
  # Let's avoid those ARM64 runners for now
@ -56,7 +57,7 @@ docker-run:
    parallel:
        matrix:
            - MATRIX_NAME: [
-                "de", "discover", "igv",
+                "de-poscounts-fallback", "discover", "igv",
                "smoke_discover", "smoke_fixed", "smoke_direct_rna", "smoke_de",
                "no_annotation", "invalid_mode", "conflicting_flags"
            ]
@ -70,18 +71,17 @@ docker-run:
          variables:
              NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
              NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq \
-                  --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf --cdna_preprocess --pychopper_backend phmm"
+                  --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf"
-              NF_IGNORE_PROCESSES: preprocess_reads,filter_unstranded_annotation,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
+              NF_IGNORE_PROCESSES: filter_unstranded_annotation,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
        - if: $MATRIX_NAME == "no_annotation"
          variables:
              NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf  ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
              NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq \
                  --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa \
                  -c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
              NF_IGNORE_PROCESSES: preprocess_reads
              ASSERT_NEXTFLOW_FAILURE: "1"
              ASSERT_NEXTFLOW_FAILURE_REXP: "Missing required parameter: --ref_annotation"
-        - if: $MATRIX_NAME == "de"
+        - if: $MATRIX_NAME == "de-poscounts-fallback"
          variables:
              NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
              NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
@ -90,7 +90,9 @@ docker-run:
                --ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gtf \
                --direct_rna --minimap2_index_opts '-k 15'  --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression/sample_sheet.csv \
                -c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
-              NF_IGNORE_PROCESSES: preprocess_reads,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
+              NF_IGNORE_PROCESSES: faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
              AFTER_NEXTFLOW_CMD: >
                  grep -Eq '"deseq2_size_factor_method": "poscounts"' ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json
        - if: $MATRIX_NAME == "only_differential_expression"
          variables:
              NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/  &&  wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
@ -103,13 +105,12 @@ docker-run:
                --igv \
                -c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
              NF_IGNORE_PROCESSES: >
-                preprocess_reads,gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,decompress_ref,
+                gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,decompress_ref,
                build_minimap_index,validate_ref_annotation,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
        - if: $MATRIX_NAME == "smoke_discover"
          variables:
              NF_BEFORE_SCRIPT: ":"
              NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf"
              NF_IGNORE_PROCESSES: preprocess_reads
              AFTER_NEXTFLOW_CMD: >
                  test -f ${CI_PROJECT_NAME}/cohort/transcripts.gtf &&
                  test -f ${CI_PROJECT_NAME}/cohort/cohort.transcriptome.fa &&
@ -118,7 +119,6 @@ docker-run:
          variables:
              NF_BEFORE_SCRIPT: ":"
              NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --transcriptome_mode fixed_annotation"
              NF_IGNORE_PROCESSES: preprocess_reads
              AFTER_NEXTFLOW_CMD: >
                  test -f ${CI_PROJECT_NAME}/cohort/transcripts.gtf &&
                  test -f ${CI_PROJECT_NAME}/cohort/transcript_counts.tsv
@ -126,15 +126,13 @@ docker-run:
          variables:
              NF_BEFORE_SCRIPT: ":"
              NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --direct_rna"
              NF_IGNORE_PROCESSES: preprocess_reads
              AFTER_NEXTFLOW_CMD: >
-                  test -f ${CI_PROJECT_NAME}/cohort/alignments/sampleA.aligned.sorted.bam &&
+                  test -f ${CI_PROJECT_NAME}/cohort/alignments/sampleA/reads.bam &&
                  test -f ${CI_PROJECT_NAME}/cohort/sqanti_cohort/classification_summary.tsv
        - if: $MATRIX_NAME == "smoke_de"
          variables:
              NF_BEFORE_SCRIPT: ":"
              NF_WORKFLOW_OPTS: "--fastq test_data/smoke/de --sample_sheet test_data/smoke/sample_sheet_de.csv --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --de_analysis --reference_level control --covariates batch"
              NF_IGNORE_PROCESSES: preprocess_reads
              AFTER_NEXTFLOW_CMD: >
                  test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dge.tsv &&
                  test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dtu_transcript.tsv
@ -142,12 +140,5 @@ docker-run:
          variables:
              NF_BEFORE_SCRIPT: ":"
              NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --transcriptome_mode nonsense"
              NF_IGNORE_PROCESSES: preprocess_reads
              ASSERT_NEXTFLOW_FAILURE: "1"
              ASSERT_NEXTFLOW_FAILURE_REXP: "nonsense is not a valid choice"
        - if: $MATRIX_NAME == "conflicting_flags"
          variables:
              NF_BEFORE_SCRIPT: ":"
              NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --direct_rna --cdna_preprocess"
              ASSERT_NEXTFLOW_FAILURE: "1"
              ASSERT_NEXTFLOW_FAILURE_REXP: ".*--cdna_preprocess cannot be used together with --direct_rna"
--- a/README.md
+++ b/README.md
@ -18,7 +18,6 @@ analysis.
 The workflow supports:
 + transcript identification from either cDNA or direct RNA reads
 + optional cDNA preprocessing with `pychopper` before alignment
 + transcript discovery guided by a supplied genome and annotation
 + quantification against a supplied reference annotation
 + optional transcript classification and QC with `SQANTI3`
@ -211,15 +210,7 @@ barcode04,treated_rep2,test_sample,treated,b2
 This example is suitable for a multiplexed run and also satisfies the minimum
 requirements for a two-group DE/DTU comparison.
-### 3. Optional cDNA preprocessing
+### 3. Genome alignment
 When `--cdna_preprocess` is enabled for cDNA libraries, the workflow runs
 `pychopper` before alignment to classify, orient, and trim full-length reads.
 This preprocessing stage is controlled by `--cdna_kit`,
 `--pychopper_backend`, and optional extra `--pychopper_opts`, and its outputs
 are published alongside the ingress results for each sample.
 ### 4. Genome alignment
 Each sample is aligned to the supplied reference genome with
 [`minimap2`](https://github.com/lh3/minimap2), then sorted and indexed with
@ -227,14 +218,14 @@ Each sample is aligned to the supplied reference genome with
 `cohort/alignments/` are the main alignment files used for transcriptome
 analysis, optional `SQANTI3` QC, and optional IGV viewing.
-### 5. Cohort transcriptome construction
+### 4. Cohort transcriptome construction
 All aligned samples are analysed together with `bambu` to produce the primary
 cohort transcriptome, transcript counts, gene counts, and the `RDS` objects used
 for downstream differential analysis. This shared model is the main cohort-level
 result and is published under `cohort/`.
-### 6. Independent per-sample transcriptomes
+### 5. Independent per-sample transcriptomes
 Each sample is also processed separately with `bambu` so the workflow produces
 sample-specific GTF, FASTA, count tables, and metadata under
@ -242,7 +233,7 @@ sample-specific GTF, FASTA, count tables, and metadata under
 specific transcript models without changing the shared cohort transcriptome used
 for DE/DTU.
-### 7. Transcript sequence generation and QC
+### 6. Transcript sequence generation and QC
 Transcript FASTA files are derived from GTF plus genome using `gffread`.
 When `--skip_sqanti` is not set, `SQANTI3` classifies the cohort and per-sample
@ -250,7 +241,7 @@ transcriptomes and produces structural QC summaries. The cohort `SQANTI3`
 results live under `cohort/sqanti_cohort/`, while per-sample `SQANTI3`
 directories are published under `samples/<alias>/<alias>_sqanti/`.
-### 8. Optional DE and DTU analysis
+### 7. Optional DE and DTU analysis
 When `--de_analysis` is enabled, the workflow checks the experimental design,
 runs `DESeq2` for differential gene expression, and runs `DEXSeq` for
@ -258,7 +249,7 @@ differential transcript usage. These analyses use the shared `bambu` outputs
 and the design columns in the sample sheet, and each comparison is written to
 its own subdirectory under `de_analysis/<contrast>/`.
-### 9. What you need to provide
+### 8. What you need to provide
 The workflow's analysis is controlled by a user provided genome, annotation, and
 `bambu` mode.
@ -270,13 +261,10 @@ The workflow's analysis is controlled by a user provided genome, annotation, and
 * `--ref_transcriptome` has been removed; if you want annotation-based
  quantification, use `--transcriptome_mode fixed_annotation` together with
  `--ref_genome` and `--ref_annotation`
 * `--cdna_preprocess` enables the `pychopper` cDNA preprocessing stage, and
  `--cdna_kit`, `--pychopper_backend`, and `--pychopper_opts` control that stage
 * when `--de_analysis` is enabled, the sample sheet must contain `alias`, the
  primary condition column, and any requested columns named in `--covariates`
 * `--cdna_preprocess` must not be combined with `--direct_rna`
-### 10. How to read the output folder
+### 9. How to read the output folder
 The published outputs are organised around a small number of top-level
 directories:
@ -315,7 +303,6 @@ directories:
 | ref_annotation | string | Reference transcript annotation in GTF or GFF format. | Required in both discover and fixed_annotation modes. |  |
 | transcriptome_mode | string | How bambu should prepare the transcriptome model. | Use discover for reference-guided transcript discovery and quantification, or fixed_annotation for quantification only against the supplied annotation. | discover |
 | direct_rna | boolean | Set this for direct RNA sequencing libraries. |  | False |
 | cdna_preprocess | boolean | Apply the optional cDNA preprocessing stage before alignment. | This is only relevant for cDNA libraries and must not be combined with `--direct_rna`. | False |
 ### Sample Options
@ -351,9 +338,6 @@ directories:
 | threads | integer | Thread count to use for the core workflow processes. |  | 4 |
 | minimap2_opts | string | Extra command-line options to pass to minimap2. |  |  |
 | ndr | number | Optional bambu novel discovery rate override. |  |  |
 | cdna_kit | string | ONT cDNA kit identifier used for pychopper preprocessing. | Only used when `--cdna_preprocess` is enabled. The workflow derives the pychopper kit code from this value. | SQK-PCS109 |
 | pychopper_backend | string | Primer-detection backend to use for pychopper preprocessing. |  | edlib |
 | pychopper_opts | string | Extra command-line options to pass to pychopper. |  |  |
 | skip_sqanti | boolean | Skip SQANTI3 transcript classification and QC. |  | False |
 | sqanti_skip_orf | boolean | Skip ORF prediction during SQANTI3 QC. |  | True |
 | sqanti_extra_args | string | Extra command-line options to pass to SQANTI3. |  |  |
@ -374,9 +358,9 @@ Output files may be aggregated including information for all samples or provided
 | Per-read stats | ingress_results/{{ alias }}/fastcat_stats/per-read-stats.tsv.gz | Read statistics for individual reads in a sample, when this output is enabled. | per-sample |
 | Ingress reads | ingress_results/{{ alias }}/seqs.fastq.gz | Reads prepared from the input data for downstream analysis. | per-sample |
 | Ingress metadata | ingress_results/{{ alias }}/metamap.json | Per-sample metadata used by the workflow. | per-sample |
-| Aligned BAM | cohort/alignments/{{ alias }}.aligned.sorted.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
+| Aligned BAM | cohort/alignments/{{ alias }}/reads.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
-| Aligned BAM index | cohort/alignments/{{ alias }}.aligned.sorted.bam.bai | Index for the aligned BAM. | per-sample |
+| Aligned BAM index | cohort/alignments/{{ alias }}/reads.bam.bai | Index for the aligned BAM. | per-sample |
-| Alignment summary | cohort/alignments/{{ alias }}.flagstat.txt | samtools flagstat output for the aligned BAM. | per-sample |
+| Alignment summary | cohort/alignments/{{ alias }}/bamstats.flagstat.tsv | bamstats flagstat summary for the aligned BAM. | per-sample |
 | Reference and annotation preparation summary | cohort/reference/annotation_reference_summary.json | Summary of reference and annotation preparation, including seqname overlap, build/provider hints, and excluded unstranded annotation counts. | aggregated |
 | Excluded unstranded annotation records | cohort/reference/unstranded_annotation.gtf | Full set of annotation records excluded because their strand was not '+' or '-'. Present only when unstranded records are found. | aggregated |
 | Cohort transcriptome GTF | cohort/transcripts.gtf | Joint bambu transcript model used as the primary cohort transcriptome. | aggregated |
@ -429,7 +413,6 @@ Find related RNA and cDNA sequencing protocols in the
  primary condition column, and any columns named in `--covariates`.
 + DE/DTU requires at least two condition levels and at least two samples per
  level.
 + `--cdna_preprocess` must not be combined with `--direct_rna`.
 + See how to interpret common Nextflow exit codes
  [here](https://labs.epi2me.io/trouble-shooting/).
@ -464,17 +447,6 @@ setting that chose how the workflow behaved. In the current version,
 The options `--transcriptome_mode discover` or `--transcriptome_mode
 fixed_annotation` should be used to choose between the modes of operation.
 #### I expected pychopper-style preprocessing for cDNA libraries
 cDNA preprocessing is available through
 `--cdna_preprocess`, which runs `pychopper` before alignment. The
 parameters for that stage are `--cdna_kit`, `--pychopper_backend`, and
 `--pychopper_opts`.
 Enable `--cdna_preprocess` for cDNA libraries when you want
 `pychopper` preprocessing, or use `--direct_rna` without
 `--cdna_preprocess` for direct RNA data.
 #### I cannot find the old flat DE output files
 In the previous workflow version, DE files appeared directly under `de_analysis/`.
@ -491,9 +463,7 @@ In the current versions, the output folder is organised around `ingress_results/
 `igv_reference/`.
 Primary shared transcriptome results are under `cohort/`, and
-`samples/<alias>/` for sample-specific models. If `--cdna_preprocess` is enabled,
+`samples/<alias>/` for sample-specific models.
 look under `ingress_results/<alias>/` for the `pychopper` outputs associated
 with each sample.
 #### My DE/DTU run fails because of `--sample_sheet`
@ -507,11 +477,6 @@ are required, and each level must contain at least two samples.
 What to change: verify the sample sheet columns first, then check
 `--condition_column`, `--covariates`, and `--reference_level`.
 #### I combined `--direct_rna` with `--cdna_preprocess`
 This is not a valid combination of parameters. Use `--direct_rna` for direct RNA
 libraries and leave `--cdna_preprocess` unset.
 #### I hit genome/annotation validation or strand-related annotation warnings
 The workflow validates that the annotation and genome share
--- a/bin/workflow_glue/models/common.py
+++ b/bin/workflow_glue/models/common.py
@ -206,10 +206,11 @@ class Sample:
        self.additional_identifiers.append(sample_identifier)
        return self.additional_identifiers
    # TODO: Integrate this into WorkflowBaseModel
    def to_json(self, filename):
        """Save class as JSON."""
        with open(filename, 'w') as f:
-            json.dump(asdict(self), f, default=str, indent=2)
+            json.dump(asdict(self), f, default=str, indent=2, allow_nan=False)
    def get_reportable_qc_status(self, max_criteria=4):
        """Store global status of the sample and list of QC criteria to show.
@ -364,7 +365,8 @@ class WorkflowResult(WorkflowBaseModel):
        self.versions = versions
        return self.versions
    # TODO: Integrate this into WorkflowBaseModel
    def to_json(self, filename):
        """Save class as JSON."""
        with open(filename, 'w') as f:
-            json.dump(asdict(self), f, default=str, indent=2)
+            json.dump(asdict(self), f, default=str, indent=2, allow_nan=False)
--- a/bin/workflow_glue/report.py
+++ b/bin/workflow_glue/report.py
@ -81,17 +81,6 @@ def _sqanti_tables(sqanti_dir):
    return tables
 def _pychopper_tables(pychopper_dir):
    tables = {}
    for summary in sorted(Path(pychopper_dir).rglob("pychopper_summary.tsv")):
        table = _read_table(summary)
        if table is None or table.empty:
            continue
        label = summary.parent.name.replace("_pychopper_output", "")
        tables[label] = table
    return tables
 def _top_results(de_dir, filename, n=20):
    tables = {}
    for contrast_dir in sorted(Path(de_dir).iterdir()):
@ -193,6 +182,7 @@ def _collect_de_method_rows(de_qc):
        fallback = contrast_data.get("deseq2_dispersion_fallback") or {}
        fallback_applied = bool(fallback.get("applied", False))
        deseq2_method = fallback.get("method_used")
        deseq2_size_factors = contrast_data.get("deseq2_size_factor_method") or "ratio"
        if not deseq2_method:
            deseq2_method = "gene-wise" if fallback_applied else "parametric"
@ -201,6 +191,7 @@ def _collect_de_method_rows(de_qc):
            deseq2_gene_wise.append(contrast_name)
        dexseq_method = contrast_data.get("dexseq_dispersion_method") or "parametric"
        dexseq_size_factors = contrast_data.get("dexseq_size_factor_method") or "ratio"
        if dexseq_method == "gene-wise":
            dexseq_gene_wise.append(contrast_name)
@ -213,11 +204,13 @@ def _collect_de_method_rows(de_qc):
        rows.append(
            {
                "Contrast": contrast_name,
                "DESeq2 size factors": deseq2_size_factors,
                "DESeq2 dispersion": (
                    f"{deseq2_method} (fallback)"
                    if fallback_applied
                    else deseq2_method
                ),
                "DEXSeq size factors": dexseq_size_factors,
                "DEXSeq dispersion": dexseq_method,
                "DEXSeq covariates dropped": (
                    ", ".join(dropped_covariates) if dropped_covariates else "none"
@ -257,10 +250,19 @@ def main(args):
            sample_names = tuple(
                item["alias"] for item in metadata if item.get("has_stats")
            )
            flagstats = tuple(
                Path(stats_dir) / "bamstats.flagstat.tsv" for stats_dir in stats
            )
            if len(stats) == 1:
                stats = stats[0]
                flagstats = flagstats[0]
                sample_names = sample_names[0] if sample_names else None
-            fastcat.SeqSummary(stats, sample_names=sample_names)
+            fastcat.SeqSummary(
                stats,
                flagstat=flagstats,
                sample_names=sample_names,
                alignment_stats=True,
            )
    with report.add_section("Sample metadata", "Samples"):
        tabs = Tabs()
@ -499,16 +501,6 @@ def main(args):
                with tabs.add_tab(stats_file.stem.replace(".flagstat", "")):
                    pre(stats_file.read_text())
    pychopper_tables = {}
    if args.pychopper_dir and Path(args.pychopper_dir).exists():
        pychopper_tables = _pychopper_tables(args.pychopper_dir)
    if pychopper_tables:
        with report.add_section("Pychopper preprocessing", "Pychopper"):
            tabs = Tabs()
            for label, table in pychopper_tables.items():
                with tabs.add_tab(label):
                    DataTable.from_pandas(table, use_index=False)
    sqanti_tables = _sqanti_tables(args.sqanti_dir)
    if sqanti_tables:
        with report.add_section("SQANTI3 classification", "SQANTI3"):
@ -853,11 +845,6 @@ def argparser():
        required=True,
        help="Per-sample output directory.",
    )
    parser.add_argument(
        "--pychopper_dir",
        default=None,
        help="Pychopper output directory.",
    )
    parser.add_argument(
        "--sqanti_dir",
        required=True,
--- a/bin/workflow_glue/tests/test_report.py
+++ b/bin/workflow_glue/tests/test_report.py
@ -217,25 +217,6 @@ def test_report_main_accepts_optional_file_sentinels(monkeypatch, tmp_path):
    assert any("GRCh38" in table.to_string() for table in tables)
 def test_pychopper_tables_uses_sample_directory_names(tmp_path):
    """Pychopper summaries should be keyed by sample alias."""
    pychopper_dir = tmp_path / "pychopper"
    sample_dir = pychopper_dir / "sampleA_pychopper_output"
    sample_dir.mkdir(parents=True)
    _write(
        sample_dir / "pychopper_summary.tsv",
        "Classification\tValue\nFull length\t10\nUnclassified\t2\n",
    )
    tables = report._pychopper_tables(pychopper_dir)
    assert list(tables) == ["sampleA"]
    assert list(tables["sampleA"]["Classification"]) == [
        "Full length",
        "Unclassified",
    ]
 def test_report_main_renders_statistical_methods_and_warnings(
    monkeypatch,
    tmp_path,
--- a/bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py
+++ b/bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py
@ -1,5 +1,5 @@
 """Check (u)BAM files for `@SQ` lines whether they are the same in all headers."""
-
+import os
 from pathlib import Path
 import sys
@ -8,58 +8,156 @@ import pysam
 from ..util import get_named_logger, wf_parser  # noqa: ABS101
 def extract_header_info(xam_file, has_reads, check_ref=False):
    """Extract header information from a BAM/CRAM file."""
    try:
        f = pysam.AlignmentFile(xam_file, check_sq=False)
    except (ValueError, IOError):
        return None, None, None, False
    with f:
        # Extract SQ lines, comparing only SN/LN/M5 elements
        # (see CW-4842 - ignore different SQ.UR values)
        sq_lines = [{
            "SN": sq["SN"],
            "LN": sq["LN"],
            "M5": sq.get("M5"),
        } for sq in f.header.get("SQ", [])]
        hd_lines = f.header.get("HD")
        xam_reflen = None
        if check_ref:
            xam_reflen = set(zip(f.references, f.lengths))
        # Whilst we have file open check for atleast 1 reads
        # unless reads already found
        if not has_reads:
            try:
                next(f.fetch(until_eof=True))
                has_reads = True
            except StopIteration:
                has_reads = False
    return sq_lines, hd_lines, xam_reflen, has_reads
 def compare_ref_lengths(xam_reflen, ref_reflen, logger):
    """Compare reference lengths from FASTA and BAM/CRAM SQ lines."""
    diff = ref_reflen.symmetric_difference(xam_reflen)
    if len(diff) > 0:
        rows = ["sequence_name sequence_length in_ref in_xam"]
        rows += [
            f"{reflen[0]} {reflen[1]} "
            f"{'1' if reflen in ref_reflen else '0'} "
            f"{'1' if reflen in xam_reflen else '0'}"
            for reflen in diff
        ]
        logger.info("\n".join(rows) + "\n")
        logger.info(
            " [WARN] In the input XAM there is at least one (name, length)"
            " pair that does not map 1:1 between the input alignment"
            " and reference. XAM will be realigned.\n"
        )
        return False
    else:
        logger.info(
            "[OKAY] Input alignment and reference sequences map 1:1\n"
        )
        return True
 def main(args):
    """Run the entry point."""
    logger = get_named_logger("checkBamHdr")
    # Check the reference if given.
    ref_reflen = None
    if args.ref:
        try:
            ref = pysam.FastaFile(
                    args.ref,
                    filepath_index=args.ref_idx)
        except (ValueError, IOError):
            logger.info(
                f"[FAIL] {args.ref} reference input file could not be"
                " read. Is it in the right format?"
            )
            sys.exit(os.EX_NOINPUT)
        # get ref length for comparison with BAM/CRAM SQ lines
        ref_reflen = set(zip(ref.references, ref.lengths))
    if not args.input_path.is_dir():
        raise ValueError(f"Input path '{args.input_path}' must be a directory.")
-    target_files = list(args.input_path.glob("*"))
+    target_files = sorted(list(args.input_path.glob("*")))
    if not target_files:
        raise ValueError(f"No files found in input directory '{args.input_path}'.")
    # Loop over target files and check if there are `@SQ` lines in all headers or not.
-    # Set `is_unaligned` accordingly. If there are mixed headers (either with some files
+    # Set `is_unaligned` accordingly.
-    # containing `@SQ` lines and some not or with different files containing different
+    # Detect mixed headers by comparing @SQ lines between files.
-    # `@SQ` lines), set `mixed_headers` to `True`.
+    # Set mixed_sq_headers to True if they differ across files
    # (either some files lack @SQ lines, or files have different @SQ lines).
    # Also check if there is the SO line, to validate whether the file is (un)sorted.
    first_sq_lines = None
-    mixed_headers = False
+    mixed_sq_headers = False
    sorted_xam = False
    xam_reflen = None
    requires_realign = False
    any_sq_lines = False  # Track if any file has @SQ lines
    has_reads = False
    for xam_file in target_files:
-        # get the `@SQ` and `@HD` lines in the header
+        sq_lines, hd_lines, xam_reflen, has_reads = extract_header_info(
-        with pysam.AlignmentFile(xam_file, check_sq=False) as f:
+            xam_file, has_reads, check_ref=bool(args.ref))
-            # compare only the SN/LN/M5 elements of SQ to avoid labelling XAM with
+
-            # same reference but different SQ.UR as mixed_header (see CW-4842)
+        if sq_lines is None:
-            sq_lines = [{
+            # File couldn't be opened
-                "SN": sq["SN"],
+            logger.error(f"Failed to open {xam_file}")
-                "LN": sq["LN"],
+            continue
-                "M5": sq.get("M5"),
+
-            } for sq in f.header.get("SQ", [])]
+        if sq_lines:
-            hd_lines = f.header.get("HD")
+            any_sq_lines = True
        # Check if it is sorted.
        # When there is more than one BAM, merging/sorting
        # will happen regardless of this flag.
        if hd_lines is not None and hd_lines.get('SO') == 'coordinate':
            sorted_xam = True
        if first_sq_lines is None:
            # this is the first file
            first_sq_lines = sq_lines
        else:
            # this is a subsequent file; check with the first `@SQ` lines
            if sq_lines != first_sq_lines:
-                mixed_headers = True
+                mixed_sq_headers = True
                break
-    # we set `is_unaligned` to `True` if there were no mixed headers and the last file
+        if ref_reflen and not mixed_sq_headers and not requires_realign:
-    # didn't have `@SQ` lines (as we can then be sure that none of the files did)
+            if sq_lines:
-    is_unaligned = not mixed_headers and not sq_lines
+                # if a reference was given and there are sq lines (indicating aligned)
-    # write `is_unaligned` and `mixed_headers` out so that they can be set as env.
+                # compare ref and xam header lengths to check alignment is to this ref
-    # variables
+                requires_realign = not compare_ref_lengths(
                    xam_reflen, ref_reflen, logger)
            else:
                requires_realign = True
    # Assumption: presence of @SQ lines indicates aligned data.
    # Set `is_unaligned` to `True` if there were no mixed headers and none of the
    # files had `@SQ` lines, or if found to be incorrectly aligned based on
    # SQ/ref length comparison.
    is_unaligned = (not mixed_sq_headers and not any_sq_lines) or requires_realign
    # write `is_unaligned`, `is_sorted` and `mixed_sq_headers`
    # out so that they can be set as env.
    # variables and handled downstream.
    sys.stdout.write(
        f"IS_UNALIGNED={int(is_unaligned)};" +
-        f"MIXED_HEADERS={int(mixed_headers)};" +
+        f"MIXED_SQ_HEADERS={int(mixed_sq_headers)};" +
-        f"IS_SORTED={int(sorted_xam)}"
+        f"IS_SORTED={int(sorted_xam)};" +
        f"HAS_READS={int(has_reads)};"
    )
    logger.info(f"Checked (u)BAM headers in '{args.input_path}'.")
@ -68,4 +166,9 @@ def argparser():
    """Argument parser for entrypoint."""
    parser = wf_parser("check_bam_headers_in_dir")
    parser.add_argument("input_path", type=Path, help="Path to target directory")
    parser.add_argument(
        "--ref", type=Path, help="Optional reference FASTA file", required=False)
    parser.add_argument(
        "--ref_idx", type=Path,
        help="Optional reference index FASTA file required for pysam", required=False)
    return parser
--- a/bin/workflow_glue_r/R/bambu.R
+++ b/bin/workflow_glue_r/R/bambu.R
@ -1,8 +1,7 @@
 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, "--bams", help = "Comma-separated BAM paths.")
-    parser <- argparser::add_argument(parser, "--bam_path", help = "Path to a single BAM file.")
+    parser <- argparser::add_argument(parser, "--aliases", help = "Comma-separated aliases for --bams.")
    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.")
@ -31,8 +30,11 @@ bambu_arg_parser <- function() {
 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)) {
+    if (workflow_glue_r_arg_missing(argv$bams)) {
-        stop("Provide exactly one of --bam_dir or --bam_path.", call. = FALSE)
+        stop("Missing required arguments: --bams", call. = FALSE)
    }
    if (workflow_glue_r_arg_missing(argv$aliases)) {
        stop("Missing required arguments: --aliases", call. = FALSE)
    }
    if (!argv$transcriptome_mode %in% c("discover", "fixed_annotation")) {
@ -52,16 +54,9 @@ bambu_validate_args <- function(argv) {
    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,
+    bamfile_list_ctor = Rsamtools::BamFileList
    list_files_fn = base::list.files
 ) {
    sample_df <- NULL
    if (!workflow_glue_r_arg_missing(argv$sample_sheet)) {
@ -81,19 +76,14 @@ bambu_resolve_inputs <- function(
        }
    }
-    if (!workflow_glue_r_arg_missing(argv$bam_dir)) {
+    bam_paths <- workflow_glue_r_parse_csv_list(argv$bams)
-        bam_paths <- sort(list_files_fn(argv$bam_dir, pattern = "\\.bam$", full.names = TRUE))
+    aliases <- workflow_glue_r_parse_csv_list(argv$aliases)
-        if (length(bam_paths) < 1) {
+
-            stop("No BAM files were found in bam_dir.", call. = FALSE)
+    if (length(bam_paths) < 1) {
-        }
+        stop("No BAM files were provided in --bams.", call. = FALSE)
-        aliases <- unname(vapply(bam_paths, bambu_strip_alias, character(1)))
+    }
-    } else {
+    if (length(aliases) != length(bam_paths)) {
-        bam_paths <- argv$bam_path
+        stop("Provide one alias per BAM in --bams.", call. = FALSE)
        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)])
@ -435,8 +425,7 @@ main_run_bambu <- function(
    prepare_annotations_fn = bambu::prepareAnnotations,
    gene_expression_fn = bambu::transcriptToGeneExpression,
    write_gtf_fn = bambu::writeToGTF,
-    bamfile_list_ctor = Rsamtools::BamFileList,
+    bamfile_list_ctor = Rsamtools::BamFileList
    list_files_fn = base::list.files
 ) {
    set.seed(42)
    suppressPackageStartupMessages({
@ -449,8 +438,7 @@ main_run_bambu <- function(
    inputs <- bambu_resolve_inputs(
        argv,
-        bamfile_list_ctor = bamfile_list_ctor,
+        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)
--- a/bin/workflow_glue_r/R/de_analysis.R
+++ b/bin/workflow_glue_r/R/de_analysis.R
@ -149,20 +149,38 @@ de_extract_disp_gene_est <- function(object) {
    S4Vectors::mcols(object)$dispGeneEst
 }
 # DESeq2's default geometric-mean size-factor estimator is undefined when
 # every gene has at least one zero across samples, so use poscounts then.
 de_choose_size_factor_type <- function(count_mat, context_label = "Count matrix") {
    if (all(rowSums(count_mat == 0) > 0)) {
        warning(
            context_label,
            " has every gene containing at least one zero; using DESeq2 size-factor estimation with sfType='poscounts'."
        )
        return("poscounts")
    }
    "ratio"
 }
 de_run_deseq_with_fallback <- function(
    dds,
    contrast_name,
    out_dir
 ) {
    sf_type <- de_choose_size_factor_type(
        DESeq2::counts(dds),
        context_label = paste0("DGE count matrix for ", contrast_name)
    )
    fallback_info <- list(
        applied = FALSE,
        method_used = "parametric",
        reason = NULL,
-        diagnostic_file = NULL
+        diagnostic_file = NULL,
        size_factor_type = sf_type
    )
    de_out <- tryCatch(
-        DESeq2::DESeq(dds, quiet = TRUE),
+        DESeq2::DESeq(dds, quiet = TRUE, sfType = sf_type),
        error = function(err) {
            if (!grepl(
                "all gene-wise dispersion estimates are within 2 orders of magnitude",
@ -184,7 +202,7 @@ de_run_deseq_with_fallback <- function(
                "Results will have reduced power and wider confidence intervals."
            )
-            dds <- DESeq2::estimateSizeFactors(dds)
+            dds <- DESeq2::estimateSizeFactors(dds, type = sf_type)
            dds <- DESeq2::estimateDispersionsGeneEst(dds)
            dds <- de_set_dispersions(dds, de_extract_disp_gene_est(dds))
@ -242,7 +260,8 @@ de_run_deseq_with_fallback <- function(
                applied = TRUE,
                method_used = "gene-wise",
                reason = conditionMessage(err),
-                diagnostic_file = basename(diag_file)
+                diagnostic_file = basename(diag_file),
                size_factor_type = sf_type
            )
            DESeq2::nbinomWaldTest(dds)
@ -391,6 +410,17 @@ de_run_deseq2_result <- function(
    )
 }
 de_estimate_size_factors_for_dexseq <- function(dxd, count_mat, coldata, sf_type) {
    sf_dds <- DESeq2::DESeqDataSetFromMatrix(
        countData = round(count_mat),
        colData = coldata,
        design = ~ 1
    )
    sf_dds <- DESeq2::estimateSizeFactors(sf_dds, type = sf_type)
    DESeq2::sizeFactors(dxd) <- DESeq2::sizeFactors(sf_dds)
    dxd
 }
 de_run_dexseq_result <- function(
    tx_counts,
    tx_meta,
@ -424,7 +454,16 @@ de_run_dexseq_result <- function(
                featureID = tx_meta$TXNAME,
                groupID = tx_meta$GENEID
            )
-            dxd <- DESeq2::estimateSizeFactors(dxd)
+            dexseq_sf_type <- de_choose_size_factor_type(
                round(tx_counts),
                context_label = "DEXSeq transcript count matrix"
            )
            dxd <- de_estimate_size_factors_for_dexseq(
                dxd,
                tx_counts,
                coldata,
                dexseq_sf_type
            )
            dispersion_result <- de_estimate_dispersions_with_fallback(
                dxd,
                "DEXSeq",
@ -446,7 +485,8 @@ de_run_dexseq_result <- function(
                dxd = dxd,
                dxr = dxr,
                dexseq_dispersion_method = dispersion_method,
-                dexseq_dispersion_reason = dispersion_reason
+                dexseq_dispersion_reason = dispersion_reason,
                dexseq_size_factor_type = dexseq_sf_type
            )
        }, error = function(err) {
            if (length(active_covariates) == 0 || !grepl(
@ -592,12 +632,14 @@ main_run_de_analysis <- function(argv) {
            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),
            deseq2_size_factor_method = "ratio",
            deseq2_dispersion_fallback = list(
                applied = FALSE,
                method_used = "parametric",
                reason = NULL,
                diagnostic_file = NULL
            ),
            dexseq_size_factor_method = "ratio",
            dexseq_dispersion_method = "parametric",
            dexseq_covariates_dropped = list()
        )
@ -632,6 +674,9 @@ main_run_de_analysis <- function(argv) {
            if (is.null(fallback_method) || identical(fallback_method, "")) {
                fallback_method <- if (fallback_applied) "gene-wise" else "parametric"
            }
            if (!is.null(fallback$size_factor_type) && !identical(fallback$size_factor_type, "")) {
                contrast_qc$deseq2_size_factor_method <- fallback$size_factor_type
            }
            contrast_qc$deseq2_dispersion_fallback <- list(
                applied = fallback_applied,
                method_used = fallback_method,
@ -737,6 +782,9 @@ main_run_de_analysis <- function(argv) {
            contrast_qc$dtu_significant_transcripts <- 0
            contrast_qc$dtu_significant_genes <- 0
        } else {
            if (!is.null(dex_res$dexseq_size_factor_type) && !identical(dex_res$dexseq_size_factor_type, "")) {
                contrast_qc$dexseq_size_factor_method <- dex_res$dexseq_size_factor_type
            }
            if (!is.null(dex_res$dexseq_dispersion_method)) {
                contrast_qc$dexseq_dispersion_method <- dex_res$dexseq_dispersion_method
            }
@ -808,6 +856,7 @@ main_run_de_analysis <- function(argv) {
            "",
            "DGE Results:",
            sprintf("  Genes tested: %d", contrast_qc$genes_tested),
            sprintf("  Size factor method: %s", contrast_qc$deseq2_size_factor_method),
            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),
@ -816,6 +865,7 @@ main_run_de_analysis <- function(argv) {
            "DTU Results:",
            sprintf("  Status: %s", contrast_qc$dtu_status),
            sprintf("  Transcripts tested: %d", contrast_qc$transcripts_tested),
            sprintf("  Size factor method: %s", contrast_qc$dexseq_size_factor_method),
            if (contrast_qc$dtu_status == "SUCCESS") {
                c(
                    sprintf(
@ -842,6 +892,10 @@ main_run_de_analysis <- function(argv) {
        function(cqc) identical(cqc$deseq2_dispersion_fallback$method_used, "gene-wise"),
        de_qc_stats$contrasts
    ))
    deseq2_poscounts <- names(Filter(
        function(cqc) identical(cqc$deseq2_size_factor_method, "poscounts"),
        de_qc_stats$contrasts
    ))
    dexseq_non_parametric <- names(Filter(
        function(cqc) {
            method <- cqc$dexseq_dispersion_method
@ -853,6 +907,10 @@ main_run_de_analysis <- function(argv) {
        function(cqc) identical(cqc$dexseq_dispersion_method, "gene-wise"),
        de_qc_stats$contrasts
    ))
    dexseq_poscounts <- names(Filter(
        function(cqc) identical(cqc$dexseq_size_factor_method, "poscounts"),
        de_qc_stats$contrasts
    ))
    dexseq_covariate_drop <- names(Filter(
        function(cqc) length(cqc$dexseq_covariates_dropped) > 0,
        de_qc_stats$contrasts
@ -867,10 +925,14 @@ main_run_de_analysis <- function(argv) {
        deseq2_dispersion_fallback_contrast_names = as.list(deseq2_dispersion_fallbacks),
        deseq2_gene_wise_contrasts = length(deseq2_gene_wise),
        deseq2_gene_wise_contrast_names = as.list(deseq2_gene_wise),
        deseq2_poscounts_contrasts = length(deseq2_poscounts),
        deseq2_poscounts_contrast_names = as.list(deseq2_poscounts),
        dexseq_non_parametric_dispersion_contrasts = length(dexseq_non_parametric),
        dexseq_non_parametric_dispersion_contrast_names = as.list(dexseq_non_parametric),
        dexseq_gene_wise_dispersion_contrasts = length(dexseq_gene_wise),
        dexseq_gene_wise_dispersion_contrast_names = as.list(dexseq_gene_wise),
        dexseq_poscounts_contrasts = length(dexseq_poscounts),
        dexseq_poscounts_contrast_names = as.list(dexseq_poscounts),
        dexseq_covariate_drop_contrasts = length(dexseq_covariate_drop),
        dexseq_covariate_drop_contrast_names = as.list(dexseq_covariate_drop),
        total_covariates_dropped = total_covariates_dropped
@ -910,7 +972,9 @@ main_run_de_analysis <- function(argv) {
                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("    DESeq2 size factors: %s", cqc$deseq2_size_factor_method),
                sprintf("    DTU status: %s", cqc$dtu_status),
                sprintf("    DEXSeq size factors: %s", cqc$dexseq_size_factor_method),
                if (cqc$dtu_status == "SUCCESS") sprintf("    DTU significant: %d genes", cqc$dtu_significant_genes) else NULL
            )
        }),
--- a/bin/workflow_glue_r/tests/testthat/test_bambu.R
+++ b/bin/workflow_glue_r/tests/testthat/test_bambu.R
@ -6,32 +6,31 @@
 #' (bin/workflow_glue/prepare_annotation_reference.py) with pytest coverage.
 #' These tests focus on bambu-specific validation and integration.
-# Workflow must accept exactly one of --bam_dir or --bam_path, not both, not neither.
+# Workflow requires explicit BAM paths and aliases.
 # Fail fast with clear error rather than passing invalid inputs to bambu.
-testthat::test_that("exactly one BAM input required", {
+testthat::test_that("BAM inputs required", {
    args <- list(
        annotation = "annotation.gtf",
        genome = "genome.fa",
        out_dir = tempfile("bambu-out-"),
-        bam_dir = NULL,
+        bams = NULL,
-        bam_path = NULL,
+        aliases = NULL,
        transcriptome_mode = "discover",
        ndr = NULL
    )
    testthat::expect_error(
        bambu_validate_args(args),
-        "Provide exactly one of --bam_dir or --bam_path"
+        "Missing required arguments: --bams"
    )
-    args$bam_dir <- tempfile("bam-dir-")
+    args$bams <- "sampleA.bam"
    args$bam_path <- "sample.bam"
    testthat::expect_error(
        bambu_validate_args(args),
-        "Provide exactly one of --bam_dir or --bam_path"
+        "Missing required arguments: --aliases"
    )
-    args$bam_path <- NULL
+    args$aliases <- "sampleA"
    testthat::expect_silent(bambu_validate_args(args))
 })
@ -42,8 +41,8 @@ testthat::test_that("invalid discovery settings rejected", {
        annotation = "annotation.gtf",
        genome = "genome.fa",
        out_dir = tempfile("bambu-out-"),
-        bam_dir = tempfile("bam-dir-"),
+        bams = "sampleA.bam",
-        bam_path = NULL,
+        aliases = "sampleA",
        transcriptome_mode = "novel",
        ndr = NULL
    )
@ -72,48 +71,26 @@ testthat::test_that("invalid discovery settings rejected", {
    testthat::expect_silent(bambu_validate_args(args))
 })
-# Derive clean sample aliases from BAM filenames by removing workflow-specific suffixes.
+# Fail fast if --bams is empty rather than passing empty input to bambu.
-testthat::test_that("BAM suffixes stripped from aliases", {
+testthat::test_that("empty BAM list rejected", {
    testthat::expect_equal(
        bambu_strip_alias("/tmp/sample.aligned.sorted.bam"),
        "sample"
    )
    testthat::expect_equal(
        bambu_strip_alias("/tmp/sample.bam"),
        "sample"
    )
 })
 # Fail fast if --bam_dir contains no BAM files rather than passing empty input to bambu.
 testthat::test_that("empty BAM directory rejected", {
    bam_dir <- tempfile("empty-bam-dir-")
    dir.create(bam_dir)
    args <- list(
-        bam_dir = bam_dir,
+        bams = "",
-        bam_path = NULL,
+        aliases = "",
        sample_alias = NULL,
        sample_sheet = NULL
    )
    testthat::expect_error(
        bambu_resolve_inputs(args, bamfile_list_ctor = function(paths, yieldSize) paths),
-        "No BAM files were found in bam_dir"
+        "No BAM files were provided in --bams"
    )
 })
 # Sample aliases must be unique and explicitly specified.
 # Derived aliases (from filenames) can collide: sample.bam vs sample.aligned.sorted.bam.
 # Sample sheets must have an 'alias' column with unique values that match BAM inputs.
 testthat::test_that("unique sample aliases required", {
    bam_dir <- tempfile("ambiguous-bam-dir-")
    dir.create(bam_dir)
    file.create(file.path(bam_dir, "sample.bam"))
    file.create(file.path(bam_dir, "sample.aligned.sorted.bam"))
    args <- list(
-        bam_dir = bam_dir,
+        bams = "sampleA.bam,sampleB.bam",
-        bam_path = NULL,
+        aliases = "sampleA,sampleA",
        sample_alias = NULL,
        sample_sheet = NULL
    )
    testthat::expect_error(
@ -121,8 +98,11 @@ testthat::test_that("unique sample aliases required", {
        "BAM aliases must be unique"
    )
-    sample_bam <- tempfile(fileext = ".bam")
+    args$aliases <- "sampleA"
-    file.create(sample_bam)
+    testthat::expect_error(
        bambu_resolve_inputs(args, bamfile_list_ctor = function(paths, yieldSize) paths),
        "Provide one alias per BAM in --bams"
    )
    missing_alias_sheet <- tempfile(fileext = ".csv")
    writeLines(
        paste(
@ -133,9 +113,8 @@ testthat::test_that("unique sample aliases required", {
        missing_alias_sheet
    )
    args <- list(
-        bam_dir = NULL,
+        bams = "sampleA.bam",
-        bam_path = sample_bam,
+        aliases = "sampleA",
        sample_alias = "sampleA",
        sample_sheet = missing_alias_sheet
    )
    testthat::expect_error(
@ -161,14 +140,9 @@ testthat::test_that("unique sample aliases required", {
 })
 # Sample sheet rows must align with BAM file order.
-# If sheet lists samples in different order than filesystem, reorder sheet to match.
+# If sheet lists samples in different order than CLI aliases, reorder sheet to match.
 # If sheet is missing aliases found in BAMs, fail.
 testthat::test_that("sample sheet reordered to match BAMs", {
    bam_dir <- tempfile("bam-dir-")
    dir.create(bam_dir)
    file.create(file.path(bam_dir, "sampleA.aligned.sorted.bam"))
    file.create(file.path(bam_dir, "sampleB.bam"))
    sample_sheet <- tempfile(fileext = ".csv")
    writeLines(
        paste(
@ -181,9 +155,8 @@ testthat::test_that("sample sheet reordered to match BAMs", {
    )
    args <- list(
-        bam_dir = bam_dir,
+        bams = "sampleA.aligned.sorted.bam,sampleB.bam",
-        bam_path = NULL,
+        aliases = "sampleA,sampleB",
        sample_alias = NULL,
        sample_sheet = sample_sheet
    )
    resolved <- bambu_resolve_inputs(
@ -239,8 +212,8 @@ testthat::test_that("transcriptome mode mapped to bambu args", {
 })
 # End-to-end unit test with mocked bambu analysis function.
-# Verifies --bam_dir input resolution, sample sheet reordering, and discovery settings.
+# Verifies --bams input resolution, sample sheet reordering, and discovery settings.
-testthat::test_that("bam_dir input with discovery mode", {
+testthat::test_that("bams input with discovery mode", {
    fixture_dir <- tempfile("bambu-discover-")
    dir.create(fixture_dir)
    bam_dir <- file.path(fixture_dir, "bams")
@ -281,9 +254,8 @@ testthat::test_that("bam_dir input with discovery mode", {
        annotation = "annotation.gtf",
        genome = "genome.fa",
        out_dir = file.path(fixture_dir, "out"),
-        bam_dir = bam_dir,
+        bams = paste(c(sample_a, sample_b), collapse = ","),
-        bam_path = NULL,
+        aliases = "sampleA,sampleB",
        sample_alias = NULL,
        sample_sheet = sample_sheet,
        transcriptome_mode = "discover",
        ndr = 0.25,
@ -492,10 +464,10 @@ testthat::test_that("QC stats match filtered results", {
 #
 # End-to-end tests with real bambu library (not mocked):
 # - Build BAMs from committed fixtures (reference.fa, annotation.gtf, reads.fastq)
-# - Run `supeRglue bambu` with both --bam_path and --bam_dir inputs
+# - Run `supeRglue bambu` with --bams input
 # - Verify output files exist and contain data for downstream workflow steps
-testthat::test_that("CLI single BAM input with fixed annotation", {
+testthat::test_that("CLI single BAM in bams input with fixed annotation", {
    fixture_dir <- tempfile("bambu-cli-")
    dir.create(fixture_dir)
@ -510,8 +482,8 @@ testthat::test_that("CLI single BAM input with fixed annotation", {
        "supeRglue",
        c(
            "bambu",
-            "--bam_path", bam_path,
+            "--bams", bam_path,
-            "--sample_alias", "sampleA",
+            "--aliases", "sampleA",
            "--sample_sheet", sample_sheet,
            "--annotation", annotation,
            "--genome", reference,
@ -546,7 +518,7 @@ testthat::test_that("CLI single BAM input with fixed annotation", {
    testthat::expect_gt(nrow(gene_counts), 0)
 })
-testthat::test_that("CLI bam_dir input preserves sample order", {
+testthat::test_that("CLI bams input preserves sample order", {
    fixture_dir <- tempfile("bambu-cli-dir-")
    dir.create(fixture_dir)
    bam_dir <- file.path(fixture_dir, "bams")
@ -573,7 +545,14 @@ testthat::test_that("CLI bam_dir input preserves sample order", {
        "supeRglue",
        c(
            "bambu",
-            "--bam_dir", bam_dir,
+            "--bams", paste(
                c(
                    file.path(bam_dir, "sampleA.aligned.sorted.bam"),
                    file.path(bam_dir, "sampleB.aligned.sorted.bam")
                ),
                collapse = ","
            ),
            "--aliases", "sampleA,sampleB",
            "--sample_sheet", sample_sheet,
            "--annotation", annotation,
            "--genome", reference,
--- a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R
+++ b/bin/workflow_glue_r/tests/testthat/test_de_analysis.R
@ -388,6 +388,24 @@ testthat::test_that("underspecified designs rejected", {
 # Fixture-driven tests for DESeq2/DEXSeq helper behaviour and metadata.
 # These avoid dependency injection and exercise the real package code paths.
 testthat::test_that("de_choose_size_factor_type uses poscounts for zero-heavy matrices", {
    count_mat <- matrix(
        c(
            0, 4, 5,
            3, 0, 2,
            7, 1, 0
        ),
        nrow = 3,
        byrow = TRUE
    )
    testthat::expect_warning(
        sf_type <- de_choose_size_factor_type(count_mat, "test matrix"),
        "using DESeq2 size-factor estimation with sfType='poscounts'"
    )
    testthat::expect_equal(sf_type, "poscounts")
 })
 testthat::test_that("de_run_deseq_with_fallback returns structured metadata", {
    testthat::skip_if_not_installed("DESeq2")
@ -721,7 +739,9 @@ testthat::test_that("CLI integration produces expected outputs", {
    testthat::expect_true("analysis_fallbacks" %in% names(de_qc))
    testthat::expect_true("contrasts" %in% names(de_qc))
    contrast_qc <- de_qc$contrasts[["condition_treated_vs_control"]]
    testthat::expect_true("deseq2_size_factor_method" %in% names(contrast_qc))
    testthat::expect_true("deseq2_dispersion_fallback" %in% names(contrast_qc))
    testthat::expect_true("dexseq_size_factor_method" %in% names(contrast_qc))
    testthat::expect_true("dexseq_dispersion_method" %in% names(contrast_qc))
    testthat::expect_true("dexseq_covariates_dropped" %in% names(contrast_qc))
 })
--- a/docs/02_introduction.md
+++ b/docs/02_introduction.md
@ -10,7 +10,6 @@ analysis.
 The workflow supports:
 + transcript identification from either cDNA or direct RNA reads
 + optional cDNA preprocessing with `pychopper` before alignment
 + transcript discovery guided by a supplied genome and annotation
 + quantification against a supplied reference annotation
 + optional transcript classification and QC with `SQANTI3`
--- a/docs/06_pipeline_overview.md
+++ b/docs/06_pipeline_overview.md
@ -60,15 +60,7 @@ barcode04,treated_rep2,test_sample,treated,b2
 This example is suitable for a multiplexed run and also satisfies the minimum
 requirements for a two-group DE/DTU comparison.
-### 3. Optional cDNA preprocessing
+### 3. Genome alignment
 When `--cdna_preprocess` is enabled for cDNA libraries, the workflow runs
 `pychopper` before alignment to classify, orient, and trim full-length reads.
 This preprocessing stage is controlled by `--cdna_kit`,
 `--pychopper_backend`, and optional extra `--pychopper_opts`, and its outputs
 are published alongside the ingress results for each sample.
 ### 4. Genome alignment
 Each sample is aligned to the supplied reference genome with
 [`minimap2`](https://github.com/lh3/minimap2), then sorted and indexed with
@ -76,14 +68,14 @@ Each sample is aligned to the supplied reference genome with
 `cohort/alignments/` are the main alignment files used for transcriptome
 analysis, optional `SQANTI3` QC, and optional IGV viewing.
-### 5. Cohort transcriptome construction
+### 4. Cohort transcriptome construction
 All aligned samples are analysed together with `bambu` to produce the primary
 cohort transcriptome, transcript counts, gene counts, and the `RDS` objects used
 for downstream differential analysis. This shared model is the main cohort-level
 result and is published under `cohort/`.
-### 6. Independent per-sample transcriptomes
+### 5. Independent per-sample transcriptomes
 Each sample is also processed separately with `bambu` so the workflow produces
 sample-specific GTF, FASTA, count tables, and metadata under
@ -91,7 +83,7 @@ sample-specific GTF, FASTA, count tables, and metadata under
 specific transcript models without changing the shared cohort transcriptome used
 for DE/DTU.
-### 7. Transcript sequence generation and QC
+### 6. Transcript sequence generation and QC
 Transcript FASTA files are derived from GTF plus genome using `gffread`.
 When `--skip_sqanti` is not set, `SQANTI3` classifies the cohort and per-sample
@ -99,7 +91,7 @@ transcriptomes and produces structural QC summaries. The cohort `SQANTI3`
 results live under `cohort/sqanti_cohort/`, while per-sample `SQANTI3`
 directories are published under `samples/<alias>/<alias>_sqanti/`.
-### 8. Optional DE and DTU analysis
+### 7. Optional DE and DTU analysis
 When `--de_analysis` is enabled, the workflow checks the experimental design,
 runs `DESeq2` for differential gene expression, and runs `DEXSeq` for
@ -107,7 +99,7 @@ differential transcript usage. These analyses use the shared `bambu` outputs
 and the design columns in the sample sheet, and each comparison is written to
 its own subdirectory under `de_analysis/<contrast>/`.
-### 9. What you need to provide
+### 8. What you need to provide
 The workflow's analysis is controlled by a user provided genome, annotation, and
 `bambu` mode.
@ -119,13 +111,10 @@ The workflow's analysis is controlled by a user provided genome, annotation, and
 * `--ref_transcriptome` has been removed; if you want annotation-based
  quantification, use `--transcriptome_mode fixed_annotation` together with
  `--ref_genome` and `--ref_annotation`
 * `--cdna_preprocess` enables the `pychopper` cDNA preprocessing stage, and
  `--cdna_kit`, `--pychopper_backend`, and `--pychopper_opts` control that stage
 * when `--de_analysis` is enabled, the sample sheet must contain `alias`, the
  primary condition column, and any requested columns named in `--covariates`
 * `--cdna_preprocess` must not be combined with `--direct_rna`
-### 10. How to read the output folder
+### 9. How to read the output folder
 The published outputs are organised around a small number of top-level
 directories:
--- a/docs/07_input_parameters.md
+++ b/docs/07_input_parameters.md
@ -17,7 +17,6 @@
 | ref_annotation | string | Reference transcript annotation in GTF or GFF format. | Required in both discover and fixed_annotation modes. |  |
 | transcriptome_mode | string | How bambu should prepare the transcriptome model. | Use discover for reference-guided transcript discovery and quantification, or fixed_annotation for quantification only against the supplied annotation. | discover |
 | direct_rna | boolean | Set this for direct RNA sequencing libraries. |  | False |
 | cdna_preprocess | boolean | Apply the optional cDNA preprocessing stage before alignment. | This is only relevant for cDNA libraries and must not be combined with `--direct_rna`. | False |
 ### Sample Options
@ -53,9 +52,6 @@
 | threads | integer | Thread count to use for the core workflow processes. |  | 4 |
 | minimap2_opts | string | Extra command-line options to pass to minimap2. |  |  |
 | ndr | number | Optional bambu novel discovery rate override. |  |  |
 | cdna_kit | string | ONT cDNA kit identifier used for pychopper preprocessing. | Only used when `--cdna_preprocess` is enabled. The workflow derives the pychopper kit code from this value. | SQK-PCS109 |
 | pychopper_backend | string | Primer-detection backend to use for pychopper preprocessing. |  | edlib |
 | pychopper_opts | string | Extra command-line options to pass to pychopper. |  |  |
 | skip_sqanti | boolean | Skip SQANTI3 transcript classification and QC. |  | False |
 | sqanti_skip_orf | boolean | Skip ORF prediction during SQANTI3 QC. |  | True |
 | sqanti_extra_args | string | Extra command-line options to pass to SQANTI3. |  |  |
--- a/docs/08_outputs.md
+++ b/docs/08_outputs.md
@ -7,9 +7,9 @@ Output files may be aggregated including information for all samples or provided
 | Per-read stats | ingress_results/{{ alias }}/fastcat_stats/per-read-stats.tsv.gz | Read statistics for individual reads in a sample, when this output is enabled. | per-sample |
 | Ingress reads | ingress_results/{{ alias }}/seqs.fastq.gz | Reads prepared from the input data for downstream analysis. | per-sample |
 | Ingress metadata | ingress_results/{{ alias }}/metamap.json | Per-sample metadata used by the workflow. | per-sample |
-| Aligned BAM | cohort/alignments/{{ alias }}.aligned.sorted.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
+| Aligned BAM | cohort/alignments/{{ alias }}/reads.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
-| Aligned BAM index | cohort/alignments/{{ alias }}.aligned.sorted.bam.bai | Index for the aligned BAM. | per-sample |
+| Aligned BAM index | cohort/alignments/{{ alias }}/reads.bam.bai | Index for the aligned BAM. | per-sample |
-| Alignment summary | cohort/alignments/{{ alias }}.flagstat.txt | samtools flagstat output for the aligned BAM. | per-sample |
+| Alignment summary | cohort/alignments/{{ alias }}/bamstats.flagstat.tsv | bamstats flagstat summary for the aligned BAM. | per-sample |
 | Reference and annotation preparation summary | cohort/reference/annotation_reference_summary.json | Summary of reference and annotation preparation, including seqname overlap, build/provider hints, and excluded unstranded annotation counts. | aggregated |
 | Excluded unstranded annotation records | cohort/reference/unstranded_annotation.gtf | Full set of annotation records excluded because their strand was not '+' or '-'. Present only when unstranded records are found. | aggregated |
 | Cohort transcriptome GTF | cohort/transcripts.gtf | Joint bambu transcript model used as the primary cohort transcriptome. | aggregated |
--- a/docs/10_troubleshooting.md
+++ b/docs/10_troubleshooting.md
@ -5,7 +5,6 @@
  primary condition column, and any columns named in `--covariates`.
 + DE/DTU requires at least two condition levels and at least two samples per
  level.
 + `--cdna_preprocess` must not be combined with `--direct_rna`.
 + See how to interpret common Nextflow exit codes
  [here](https://labs.epi2me.io/trouble-shooting/).
@ -40,17 +39,6 @@ setting that chose how the workflow behaved. In the current version,
 The options `--transcriptome_mode discover` or `--transcriptome_mode
 fixed_annotation` should be used to choose between the modes of operation.
 #### I expected pychopper-style preprocessing for cDNA libraries
 cDNA preprocessing is available through
 `--cdna_preprocess`, which runs `pychopper` before alignment. The
 parameters for that stage are `--cdna_kit`, `--pychopper_backend`, and
 `--pychopper_opts`.
 Enable `--cdna_preprocess` for cDNA libraries when you want
 `pychopper` preprocessing, or use `--direct_rna` without
 `--cdna_preprocess` for direct RNA data.
 #### I cannot find the old flat DE output files
 In the previous workflow version, DE files appeared directly under `de_analysis/`.
@ -67,9 +55,7 @@ In the current versions, the output folder is organised around `ingress_results/
 `igv_reference/`.
 Primary shared transcriptome results are under `cohort/`, and
-`samples/<alias>/` for sample-specific models. If `--cdna_preprocess` is enabled,
+`samples/<alias>/` for sample-specific models.
 look under `ingress_results/<alias>/` for the `pychopper` outputs associated
 with each sample.
 #### My DE/DTU run fails because of `--sample_sheet`
@ -83,11 +69,6 @@ are required, and each level must contain at least two samples.
 What to change: verify the sample sheet columns first, then check
 `--condition_column`, `--covariates`, and `--reference_level`.
 #### I combined `--direct_rna` with `--cdna_preprocess`
 This is not a valid combination of parameters. Use `--direct_rna` for direct RNA
 libraries and leave `--cdna_preprocess` unset.
 #### I hit genome/annotation validation or strand-related annotation warnings
 The workflow validates that the annotation and genome share
--- a/lib/common.nf
+++ b/lib/common.nf
@ -61,3 +61,79 @@ process configure_igv {
    """
 }
 process minimap2_alignment {
    label "wf_common"
    cpus { bsargs["alignment_threads"] }
    memory { bsargs["minimap2_memory"][task.attempt - 1] }
    maxRetries { bsargs["minimap2_memory"].size() - 1 }
    errorStrategy = {task.exitStatus in [137,140] ? 'retry' : 'finish'}
    input:
        tuple val(meta), path(reads, stageAs: 'reads/*'), path(reference), path(ref_index)
        tuple val(align_ext), val(index_ext) // either [bam, bai] or [cram, crai]
        val bsargs
    output:
        tuple val(meta),
              path("reads.${align_ext}"), 
              path("reads.${align_ext}.${index_ext}"),
              path("bamstats_results"),
              emit: alignment
    script:
        String reset_cmd_body = "samtools reset -x tp,cm,s1,s2,NM,MD,AS,SA,ms,nn,ts,cg,cs,dv,de,rl"
        String fastq_cmd_body = "samtools fastq -T *"
        // Default required threads is 6
        // Samtools x3 and bamstats will all be single threaded
        def minimap2_threads = Math.max(task.cpus - 4 , 1)
        def per_read_stats_arg = bsargs["per_read_stats"] ? "| bgzip > bamstats_results/bamstats.readstats.tsv.gz" : " > /dev/null"
        def bam_input_cmd = """samtools view -H --no-PG reads/\"\$(ls reads | head -n 1)\" \
            > reads.header && samtools cat reads/* \
            | ${reset_cmd_body} --no-PG - -o -  \
            | ${fastq_cmd_body} - """
        def fastq_input_cmd = "touch reads.header && fastcat --reheader reads/*"
        def bam_or_fastq_input = params.bam ? bam_input_cmd : fastq_input_cmd
        def minimap2_opts = bsargs["minimap2_opts"]
    """
    rm -rf bamstats_results
    mkdir bamstats_results
    ${bam_or_fastq_input} \
        | minimap2 -y -t ${minimap2_threads} -a ${minimap2_opts} --cap-kalloc 100m --cap-sw-mem 50m \
            ${reference} - \
        | workflow-glue reheader_samstream reads.header \
            --insert \$'@PG\\tID:reset\\tPN:samtools\\tCL:${reset_cmd_body}' \
            --insert \$'@PG\\tID:fastq\\tPN:samtools\\tCL:${fastq_cmd_body}' \
        | samtools sort -u -O BAM - \
        | tee >(samtools view --reference ${reference} \
            --write-index -o reads.${align_ext}##idx##reads.${align_ext}.${index_ext} -) \
        | bamstats -s ${meta.alias} -u \
            -f bamstats_results/bamstats.flagstat.tsv  \
            -i bamstats_results/bamstats.runids.tsv \
            -l bamstats_results/bamstats.basecallers.tsv \
            --histograms histograms - \
            ${per_read_stats_arg}
        # post-pipe bamstats tidying
        mv histograms/* bamstats_results/
        # get n_seqs from flagstats - need to sum them up
        awk 'NR==1{for (i=1; i<=NF; i++) {ix[\$i] = i}} NR>1 {c+=\$ix["total"]} END{print c}' \
            bamstats_results/bamstats.flagstat.tsv > bamstats_results/n_seqs
        # get unique run IDs (we add `-F '\\t'` as `awk` uses any stretch of whitespace
        # as field delimiter otherwise and thus ignore empty columns)
        awk -F '\\t' '
            NR==1 {for (i=1; i<=NF; i++) {ix[\$i] = i}}
            # only print run_id if present
            NR>1 && \$ix["run_id"] != "" {print \$ix["run_id"]}
        ' bamstats_results/bamstats.runids.tsv | sort | uniq > bamstats_results/run_ids
        # get unique basecall models
        awk -F '\\t' '
           NR==1 {for (i=1; i<=NF; i++) {ix[\$i] = i}}
           # only print model if present
           NR>1 && \$ix["basecaller"] != "" {print \$ix["basecaller"]}
        ' bamstats_results/bamstats.basecallers.tsv | sort | uniq > bamstats_results/basecallers ;
    """
 }
--- a/lib/ingress.nf
+++ b/lib/ingress.nf
@ -4,6 +4,10 @@ import ArgumentParser
 N_OPEN_FILES_LIMIT = 128
 OPTIONAL_FILE = file("$projectDir/data/OPTIONAL_FILE")
 include { minimap2_alignment as bam_alignment;
          minimap2_alignment as fastq_alignment; } from './common'
 /**
 * Check if a file ends with one of the target extensions.
@ -171,6 +175,14 @@ def add_number_of_reads_to_meta(ch, String input_type_format) {
 *  - "allow_multiple_basecall_models": emit data of samples that had more than one
 *     basecall model; if this is `false`, such samples will be emitted as `[meta, null,
 *     null]`
 *  - "minimap2_memory": A list of memory options to pass to minimap2. Alignment will be retried
 *     with next in list if it fails due to memory issues. (default: ["8G", "15G", "31G"])
 *  - "minimap2_opts": string with options to pass to minimap2 (default: "-x lr:hq")
 *  - "alignment_threads": number of threads to use for alignment process (default: 6)
 *  - "output_xam_fmt": alignment output format, `bam` outputs a BAM file with index (.bam, .bai),
       `cram` outputs a CRAM file with index (.cram, .crai)
 * @param aln_ref_ch: optional channel with a reference tuple (ref, ref_idx) to align against.
 *  If provided alignment will be attempted if inputs are unaligned or not already aligned to this ref.
 * @return: channel of `[Map(alias, barcode, type, ...), Path|null, Path|null]`.
 *  The first element is a map with metadata, the second is the path to the
 *  `.fastq.gz` file with the (potentially concatenated) sequences and the third is
@ -178,7 +190,7 @@ def add_number_of_reads_to_meta(ch, String input_type_format) {
 *  `null` for sample sheet entries for which no corresponding barcode directory was
 *  found. The third element is `null` if there were no reads.
 */
-def fastq_ingress(Map arguments)
+def fastq_ingress(Map arguments, aln_ref_ch = null)
 {
    // check arguments
    Map margs = parse_arguments(
@ -194,42 +206,99 @@ def fastq_ingress(Map arguments)
    def input = get_valid_inputs(margs, fq_extensions)
-    def ch_result = fastcat(input.files.mix(input.dirs), margs, "FASTQ")
+    // Alignment with reference if provided
-  
+    if (aln_ref_ch) {
    // TODO: xam_ingress mixes in a .no_files channel here. Do we need to do the same? 
-    // The above may have returned a channel with multiple fastqs if chunking
+        if (margs["output_xam_fmt"] == "bam"){
-    // is enabled. Flatten this and add a groupKey to meta information which
+            output_xam_fmt = ["bam", "bai"]
    // states the number of sibling files. This can be later used as the key
    // for .groupTuple() on a channel in order to get all results for a sample
    // We don't decorate "alias" with a count because that messes up downstream
    // serialisation.
    // Mix in the missing files from the sample sheet
    // Add in a unique key for every emission
    def ch_spread_result = ch_result
        .mix (input.missing.map { meta, files -> [meta, files, null] })
        .map { meta, files, stats ->
            // new `arity: '1..*'` would be nice here
            files = files instanceof List ? files : [files]
            def new_keys = [
                "group_key": groupKey(meta["alias"], files.size()),
                "n_fastq": files.size()]
            def grp_index = (0..<files.size()).collect()
            [meta + new_keys, files, grp_index, stats]
        }
-        .transpose(by: [1, 2])  // spread multiple fastq files into separate emissions
+        else {
-        .map { meta, files, grp_i, stats ->
+            output_xam_fmt = ["cram", "crai"]
-            def new_keys = [
+        }
-                "group_index": "${meta["alias"]}_${grp_i}"]
+        
-            [meta + new_keys, files, stats]
+        mm2_aln = fastq_alignment(
            input.files.mix(input.dirs).combine(aln_ref_ch),
            output_xam_fmt,
            margs
        )
        // Add back samples from sample sheet with missing barcode dir
        aligned_and_missing = mm2_aln.mix(
            input.missing.map { meta, files -> [meta, files, null, null] }
        )
        // add number of reads to meta
        updated_ch = add_number_of_reads_to_meta(
            aligned_and_missing 
            | map { meta, bam, bai, stats ->
                [meta, [bam, bai], stats]
            }, 
            "xam"
        )
        // add run IDs, and basecall models to meta
        ch_result = add_run_IDs_and_basecall_models_to_meta(
            updated_ch,
            margs.allow_multiple_basecall_models
        ) | map { it -> 
            def (meta, bam, bai, stats) = it.flatten()
            def xam = meta.src_xam
            def xai = meta.src_xai
            // S3 paths are remote references; only keep local file paths for downstream processing
            if (meta.src_xam){
                xam = meta.src_xam.startsWith('s3://') ? null : meta.src_xam
                xai = meta.src_xam.startsWith('s3://') ? null : meta.src_xai
            }
            [ meta + [src_xam: xam, src_xai: xai], bam, bai, stats ]
        }
-    // add number of reads, run IDs, and basecall models to meta
+        return ch_result
-    def ch_final = add_number_of_reads_to_meta(ch_spread_result, "fastq")
+    } else {
-    ch_final = add_run_IDs_and_basecall_models_to_meta(
+        // If no alignment ref is provided, just process FASTQ as usual.
-        ch_final, margs.allow_multiple_basecall_models
+        def ch_fastcat = fastcat(input.files.mix(input.dirs), margs, "FASTQ")
-    )
+        def ch_fastcat_branched = ch_fastcat.branch { meta, files, stats ->
-    return ch_final
+            def file_list = files ? (files instanceof List ? files : [files]) : []
            no_files: file_list.size() == 0
            has_files: true
        }
        ch_result = ch_fastcat_branched.has_files
        def ch_no_files = ch_fastcat_branched.no_files
            .map { meta, files, stats -> [meta, null, null] }
        // fastcat may have returned a channel with multiple fastqs if chunking
        // is enabled. Flatten this and add a groupKey to meta information which
        // states the number of sibling files. This can be later used as the key
        // for .groupTuple() on a channel in order to get all results for a sample
        // We don't decorate "alias" with a count because that messes up downstream
        // serialisation.
        // Mix in the missing files from the sample sheet
        // Add in a unique key for every emission
        ch_spread_result = ch_result
            .mix(input.missing.map { meta, files -> [meta, files, null] })
            .mix(ch_no_files)
            .map { meta, files, stats ->
                // new `arity: '1..*'` would be nice here
                files = files instanceof List ? files : [files]
                def new_keys = [
                    "group_key": groupKey(meta["alias"], files.size()),
                    "n_fastq": files.size()]
                def grp_index = (0..<files.size()).collect()
                [meta + new_keys, files, grp_index, stats]
            }
            .transpose(by: [1, 2])  // spread multiple fastq files into separate emissions
            .map { meta, files, grp_i, stats ->
                def new_keys = [
                    "group_index": "${meta["alias"]}_${grp_i}"]
                [meta + new_keys, files, stats]
            }
        // add number of reads, run IDs, and basecall models to meta
        ch_final = add_number_of_reads_to_meta(ch_spread_result, "fastq")
        ch_final = add_run_IDs_and_basecall_models_to_meta(
            ch_final, margs.allow_multiple_basecall_models
        )
        return ch_final
    }
 }
@ -251,7 +320,7 @@ def fastq_ingress(Map arguments)
 *    directories.
 *  - "keep_unaligned": boolean whether to include uBAM files
 *  - "return_fastq": boolean whether to convert to FASTQ (this will always run
- *    `fastcat`)
+ *    `fastcat`). Cannot be set to true if an alignment reference is provided.
 *  - "fastcat_extra_args": string with extra arguments to pass to `fastcat`
 *  - "required_sample_types": list of zero or more required sample types expected to be present 
 *     in the sample sheet
@ -260,7 +329,15 @@ def fastq_ingress(Map arguments)
 *  - "fastq_chunk": null or a number of reads to place into chunked FASTQ files
 *  - "allow_multiple_basecall_models": boolean. If true, emit data of samples that had more than one
 *     basecall model; if this is `false`, such samples will be emitted as `[meta, null,
- *     null]`
+ null]`
 *  - "minimap2_memory": a list of memory options to pass to minimap2, alignment will be retried
 *     with next in list if it fails due to memory issues. (default: ["8G", "15G", "31G"])
 *  - "minimap2_opts": string with options to pass to minimap2 (default: "-x lr:hq")
 *  - "alignment_threads": number of threads to use for alignment process (default: 6)
 *  - "output_xam_fmt": alignment output format, `bam` outputs an BAM file with index (.bam, .bai),
       `cram` outputs a CRAM file with index (.cram, .crai)
 * @param aln_ref_ch: optional channel with a reference tuple (ref, ref_idx) to align against.
 *  If provided alignment will be attempted if inputs are unaligned or not already aligned to this ref.
 * @return: channel of `[Map(alias, barcode, type, ...), Path|null, Path|null]`.
 *  The first element is a map with metadata, the second is the path to the
 *  `.bam` file with the (potentially merged) sequences and the third is
@ -269,7 +346,7 @@ def fastq_ingress(Map arguments)
 *  found and for samples with only uBAM files when `keep_unaligned: false`. The third
 *  element is `null` if `bamstats` was not run.
 */
-def xam_ingress(Map arguments)
+def xam_ingress(Map arguments, aln_ref_ch = null)
 {
    // check arguments
    Map margs = parse_arguments(
@ -283,13 +360,25 @@ def xam_ingress(Map arguments)
    )
    margs["fastq_chunk"] ?= 0  // cant pass null through channel
    if (margs["return_fastq"] && aln_ref_ch){
        error "`return_fastq` ingress argument cannot be true when alignment is enabled. Alignment always produces XAM-format outputs."
    }
    // we only accept BAM or uBAM for now (i.e. no SAM or CRAM)
    ArrayList xam_extensions = [".bam", ".ubam"]
    def input = get_valid_inputs(margs, xam_extensions)
    if (aln_ref_ch) {
        aln_ref = aln_ref_ch
    } else {
        aln_ref = Channel.of( 
          tuple(OPTIONAL_FILE, OPTIONAL_FILE)
        )
    }
    // check BAM headers to see if any samples are uBAM
-    ch_result = input.dirs
+    ch_check_bams = input.dirs
    | map { meta, path -> [meta, get_target_files_in_dir(path, xam_extensions, margs)] }
    | mix(input.files)
    | map{
@ -313,21 +402,65 @@ def xam_ingress(Map arguments)
        }
        [meta + [src_xam: src_xam, src_xai: src_xai], paths]
    }
    | combine(aln_ref)
    | checkBamHeaders
-    | map { meta, paths, is_unaligned_env, mixed_headers_env, is_sorted_env ->
+    | map { meta, paths, is_unaligned_env, mixed_sq_headers_env, is_sorted_env, has_reads_env ->
        // convert the env. variables from strings ('0' or '1') into bools
        boolean is_unaligned = is_unaligned_env as int as boolean
-        boolean mixed_headers = mixed_headers_env as int as boolean
+        boolean mixed_sq_headers = mixed_sq_headers_env as int as boolean
        boolean is_sorted = is_sorted_env as int as boolean
        // if no reads, no error (multisample may have some empty bams) but do not attempt alignment
        boolean has_reads = has_reads_env as int as boolean
        // throw an error if there was a sample with mixed headers
-        if (mixed_headers) {
+        if (mixed_sq_headers) {
            error "Found mixed headers in (u)BAM files of sample '${meta.alias}'."
        }
        // add `is_unaligned` to the metamap (note the use of `+` to create a copy of
        // `meta` to avoid modifying every item in the channel;
        // https://github.com/nextflow-io/nextflow/issues/2660)
-        [meta + [is_unaligned: is_unaligned, is_sorted: is_sorted], paths]
+        [meta + [is_unaligned: is_unaligned, is_sorted: is_sorted, has_reads: has_reads], paths]
    }
    // Handle alignment
    if (aln_ref_ch) {
        alignment_fork = ch_check_bams
        | branch {
            meta, paths ->
                to_align: (meta.is_unaligned == true) && (meta.has_reads == true)
                noalign: true
        }
        if (margs["output_xam_fmt"] == "bam"){
            output_xam_fmt = ["bam", "bai"]
        }
        else {
            output_xam_fmt = ["cram", "crai"]
        }
        mm2_aln = bam_alignment(
            alignment_fork.to_align.combine(aln_ref_ch),
            output_xam_fmt,
            margs
        )
        // Update meta is unaligned
        mm2_aln_final = mm2_aln.alignment.map{
            meta, xam, xai, stats ->
                // remove has reads from meta as no longer required
                def newmeta = meta.findAll { k, v -> k != 'has_reads' }
                [newmeta + [is_unaligned: false], xam, xai, stats]
        }
        // Process BAM files that do not require realignment by passing them through the standard downstream steps (merging, sorting, indexing, etc.)
        ch_result_tmp = alignment_fork.noalign.map{
            meta, paths -> 
                def newmeta = meta.findAll { k, v -> k != 'has_reads' }
                [newmeta, paths]
        }
    } else {
        // If no alignment reference provided process all BAM's as usual.
        ch_result_tmp = ch_check_bams.map{ meta, paths -> 
            def newmeta = meta.findAll { k, v -> k != 'has_reads' }
            [newmeta, paths]}
    }
    ch_result = ch_result_tmp 
    | branch { meta, paths ->
        // set `paths` to `null` for uBAM samples if unallowed (they will be added to
        // the results channel in shape of `[meta, null]` at the end of the function
@ -485,7 +618,8 @@ def xam_ingress(Map arguments)
    }
    // Combine all possible inputs
-    ch_result = ch_missing | mix(
+    ch_result = ch_missing 
    | mix(
        ch_no_op,
        ch_indexed,
        ch_merged,
@ -514,6 +648,11 @@ def xam_ingress(Map arguments)
        ch_result.is_null.map{it + [null]}
    )
    // Add back aligned if alignment reference provided
    if (aln_ref_ch) {
        ch_result = ch_result.mix(mm2_aln_final)
    }
    // Remove metadata that are unnecessary downstream:
    // meta.src_xai: not needed, as it will be part of the channel as a file
    // meta.is_sorted: if data are aligned, they will also be sorted/indexed
@ -550,13 +689,14 @@ def xam_ingress(Map arguments)
        meta, bam, bai, stats ->
        def xam = meta.src_xam
        def xai = meta.src_xai
        // S3 paths are remote references;
        // only keep local file paths for downstream processing
        if (meta.src_xam){
            xam = meta.src_xam.startsWith('s3://') ? null : meta.src_xam
            xai = meta.src_xam.startsWith('s3://') ? null : meta.src_xai
        }
        [ meta + [src_xam: xam, src_xai: xai], bam, bai, stats ]
    }
    return ch_result
 }
@ -629,18 +769,21 @@ process checkBamHeaders {
    label "wf_common"
    cpus 1
    memory "2 GB"
-    input: tuple val(meta), path("input_dir/reads*.bam")
+    input: tuple val(meta), path("input_dir/reads*.bam"), path(ref, stageAs: "ref/*"),
            path(ref_index, stageAs: "ref_index/*")
    output:
        tuple(
            val(meta),
            path("input_dir/reads*.bam", includeInputs: true),
            env(IS_UNALIGNED),
-            env(MIXED_HEADERS),
+            env(MIXED_SQ_HEADERS),
            env(IS_SORTED),
            env(HAS_READS),
        )
    script:
    String ref_arg = ref.fileName.name == OPTIONAL_FILE.name ? "" : "--ref $ref --ref_idx $ref_index"
    """
-    workflow-glue check_bam_headers_in_dir input_dir > env.vars
+    workflow-glue check_bam_headers_in_dir input_dir $ref_arg > env.vars
    source env.vars
    """
 }
@ -836,6 +979,10 @@ Map parse_arguments(String func_name, Map arguments, Map extra_kwargs=[:]) {
        "required_sample_types": [],
        "per_read_stats": false,
        "allow_multiple_basecall_models": false,
        "minimap2_memory": ["8GB", "15GB", "31GB"],
        "minimap2_opts": "-x lr:hq",
        "alignment_threads": 6,
        "output_xam_fmt": "bam" // or cram
    ]
    ArgumentParser parser = new ArgumentParser(
        args: required_args,
--- a/lib/reference.nf
+++ b/lib/reference.nf
@ -0,0 +1,246 @@
 // This subworkflow pre-processes an input reference
 // genome by preparing the appropriate indexes and cache.
 // By default, the workflow takes an input reference,
 // decompresses it if it is gzipped, and either imports
 // an existing index (if available) or generates one
 // using `samtools faidx`.
 // If requested by the user, the workflow can also:
 // 1. Generate a CRAM cache, using the `REF_PATH`
 //    environment variable. IMPORTANT: If the ingress workflow
 //    outputs CRAM files (when output_xam_fmt='cram'), a CRAM cache
 //    MUST be generated here and passed to downstream processes.
 //    Even though users cannot currently input CRAMs, the ingress
 //    workflow can generate them, and downstream processes require
 //    the cache to decompress CRAM files.
 // 2. Generate a minimap2 `.mmi` index for faster alignment.
 Map colors = NfcoreTemplate.logColours(params.monochrome_logs)
 // Argument parser
 Map parse_reference(Map arguments) {
    ArgumentParser parser = new ArgumentParser(
        args:[
        ],
        kwargs:[
            "output_cache": false,
            "output_mmi": false,
            "mmi_opts": "-x lr:hq",
            "mmi_memory": ["8GB", "15GB", "31GB"],
        ],
        name: "reference_ingress")
    return parser.parse_args(arguments)
 }
 // Process to generate the CRAM cache and
 // create the REF_PATH variable.
 // NOTE: The CRAM cache is essential for any downstream processes that
 // need to decompress or work with CRAM files. Even though users cannot
 // currently input CRAMs, if output_xam_fmt='cram' is used in ingress,
 // downstream processes will require this cache to properly handle the
 // CRAM output from ingress.
 process cram_cache {
    label "wf_common"
    cpus 1
    memory 4.GB
    input:
        path reference
    output:
        tuple path("ref_cache/"), env(REF_PATH), emit: ref_cache
    script:
    """
    # Invoke from binary installed to container PATH
    seq_cache_populate.pl -root ref_cache/ "${reference}"
    REF_PATH="ref_cache/%2s/%2s/%s"
    """
 }
 // Process to create the faidx index
 // main.nf currently handles publishing the faidx if required
 // maintain naming so it matches input ref when output
 process faidx {
    label "wf_common"
    cpus 1
    memory 4.GB
    input:
        path(ref)
    output:
        path("${ref}.fai")
    script:
    """
    samtools faidx "${ref}"
    """
 }
 // Process to create the faidx indexes for a gzipped reference
 process gz_faidx {
    publishDir "${params.out_dir}", mode: 'copy', pattern: "*"
    label "wf_common"
    cpus 1
    memory 4.GB
    // If a user provides a non-bgzipped file indexing will fail, which should be tolerated.
    // With error_strategy set to 'ignore', this will result in an empty channel 
    // that can be checked for downstream using `ifEmpty()` 
    errorStrategy 'ignore'
    input:
        path(ref)
    output:
        tuple path(ref), path("${ref}.fai"), path("${ref}.gzi")
    script:
    """
    samtools faidx "${ref}"
    """
 }
 // Decompress the reference genome
 // NOTE -f required to decompress symlink
 process decompress_ref {
    label "wf_common"
    cpus 1
    memory 4.GB
    input:
        path "ref.fa.gz"
    output:
        path "ref.fa", emit: decompressed_ref
    """
    gzip -df ref.fa.gz
    """
 }
 // Prepare minimap2 .mmi index
 process make_mmi {
    memory {  margs["mmi_memory"][task.attempt - 1] }
    maxRetries { margs["mmi_memory"].size() - 1 }
    errorStrategy 'retry'
    label "wf_common"
    cpus 3
    input:
        path("ref.fa")
        val margs
    output:
        path("ref.mmi")
    script:
    """
    minimap2 -t ${task.cpus} ${margs["mmi_opts"]} -d ref.mmi ref.fa
    """
 }
 // Workflow to prepare the reference genome and its indexes.
 /**
 * This workflow accepts either a string file path or a Nextflow channel as the reference input,
 * making it flexible for both direct file inputs and outputs from upstream processes.
 * 
 * IMPORTANT: If your workflow uses CRAM output from ingress (output_xam_fmt='cram'),
 * you must set "output_cache": true to generate the CRAM cache. This cache is required
 * by downstream processes to properly handle CRAM files, even though users cannot currently
 * input CRAMs directly.
 * 
 * inputs:
 *  - reference_path_or_channel: (String or Channel) Path to reference FASTA file or a Nextflow channel emitting a reference file
 *  - arguments: (Map) Optional configuration parameters:
 *    - "output_cache": (boolean) If true, generate a CRAM cache for the reference (default: false)
 *                      REQUIRED if using CRAM output from ingress for downstream processing.
 *    - "output_mmi": (boolean) If true, generate a minimap2 .mmi index for the reference (default: false)
 *    - "mmi_opts": (string) Options string passed to minimap2 for index creation (default: "-x lr:hq")
 *    - "mmi_memory": (list) Memory options to pass to minimap2 index; creation will be retried with next in list if it fails due to memory issues (default: ["8GB", "15GB", "31GB"])
 *
 * emits:
 *  - "ref_tuple": Tuple of (reference, fai_index)
 *  - "ref_cache": CRAM cache for the reference genome if "output_cache" is true.
 *                 Required for downstream processes when CRAM files are used.
 *  - "ref_mmi": Minimap2 index for the reference genome if "output_mmi" is true
 *  - "ref_gzidx": Gzip index for the reference genome (useful for IGV)
 */
 workflow prepare_reference {
    take:
        reference_path_or_channel
        arguments
    main:
        Map margs = parse_reference(arguments)
        // Determine if input is a string path or channel
        Boolean is_string = reference_path_or_channel in String
        // Base ref channel
        if (!is_string) {
            // It's a channel - use it directly
            ref = reference_path_or_channel
            is_compressed = false  // Assume uncompressed when coming from channel
            gzindexes = Channel.empty()
        } else {
            // It's a string path
            if (file(reference_path_or_channel).exists()){
                ref = Channel.fromPath(reference_path_or_channel)
            } else {
                throw new Exception(colors.red + "File ${reference_path_or_channel} not found." + colors.reset)
            }
            // Check if compressed
            is_compressed = reference_path_or_channel.toLowerCase().endsWith("gz")
            if (is_compressed) {
                // Define indexes names.
                String input_fai_index = "${reference_path_or_channel}.fai"
                String input_gzi_index = "${reference_path_or_channel}.gzi"
                // Decompress reference genome.
                ref = decompress_ref(ref)
                // Check whether the input gzref is indexed. If so, pass these as indexes.
                // Otherwise, generate the gzip + fai indexes for the compressed reference.
                if (file(input_fai_index).exists() && file(input_gzi_index).exists()){
                    gzindexes = Channel.fromPath(reference_path_or_channel)
                    | mix(
                        Channel.fromPath(input_fai_index),
                        Channel.fromPath(input_gzi_index)
                    )
                } else {
                    gzindexes = gz_faidx(Channel.fromPath(reference_path_or_channel))
                    gzindexes.ifEmpty{
                        if (params.containsKey("igv") && params.igv){
                            log.warn """\
                                The input reference is compressed but not with bgzip, which is required to create an index.
                                The workflow will proceed but it will not be possible to load the reference in the IGV Viewer.
                                To use the IGV Viewer, provide an uncompressed, or bgzip compressed version of the input reference next time you run the workflow.
                                """.stripIndent()
                        }
                    }
                }
            } else {
                gzindexes = Channel.empty()
            }
        }
        // Generate fai index if the file is either compressed, or if fai doesn't exists
        if (!is_compressed) {
            if (is_string && file("${reference_path_or_channel}.fai").exists()){
                ref_idx = Channel.fromPath("${reference_path_or_channel}.fai")
            } else {
                ref_idx = faidx(ref)
            }
        } else {
            ref_idx = faidx(ref)
        }
        // Combine ref and its fai index
        ref_tuple = ref.combine(ref_idx)
        // Generate CRAM cache
        if (margs.output_cache){
            cram_cache(ref)
            ref_cache = cram_cache.out.ref_cache
        } else {
            ref_cache = null
        }
        // Generate mmi index
        if (margs.output_mmi){
            ref_mmi = make_mmi(ref, margs)
        } else {
            ref_mmi = null
        }
    emit:
        ref_tuple = ref_tuple
        ref_cache = ref_cache
        ref_mmi = ref_mmi
        ref_gzidx = gzindexes
 }
--- a/main.nf
+++ b/main.nf
@ -5,6 +5,7 @@ nextflow.enable.dsl = 2
 include { fastq_ingress; xam_ingress } from './lib/ingress'
 include { getParams; configure_igv } from './lib/common'
 include { prepare_reference } from './lib/reference'
 include { transcriptome_analysis } from './subworkflows/transcriptome'
 include { differential_expression } from './subworkflows/differential_expression'
@ -29,69 +30,6 @@ process getVersions {
 }
 process collectIngressResultsInDir {
    label "wf_common"
    cpus 1
    memory "2 GB"
    input:
        tuple val(meta), path(reads, stageAs: "reads/*"), path(stats, stageAs: "stats/*")
    output:
        path "out/*", emit: out
    script:
        String outdir = "out/${meta["alias"].replaceAll("'", "'\\\\''")}"
        String meta_json = new JsonBuilder(meta).toPrettyString().replaceAll("'", "'\\\\''")
        String stats_arg = stats.fileName.name == OPTIONAL_FILE.name ? "" : stats
        def read_args = reads instanceof java.util.ArrayList ? reads.join(" ") : reads
    """
    mkdir -p '${outdir}'
    echo '${meta_json}' > metamap.json
    mv metamap.json ${read_args} ${stats_arg} '${outdir}'
    """
 }
 process preprocess_reads {
    label "wf_transcriptomes_pychopper"
    cpus { params.threads ?: 4 }
    memory "8 GB"
    input:
        tuple val(meta), path(reads, stageAs: "reads/*")
    output:
        tuple val(meta.alias), path("${meta.alias}_pychopper_output/${meta.alias}_full_length_reads.fastq"), emit: full_len_reads
        tuple val(meta.alias), path("${meta.alias}_pychopper_output"), emit: dir
    script:
        String backend = params.pychopper_backend ?: "edlib"
        String extra = params.pychopper_opts ?: ""
        String cdna_kit = params.cdna_kit ? params.cdna_kit.tokenize("-")[-1] : ""
        def read_args = reads instanceof java.util.ArrayList ? reads.join(" ") : reads
    """
    cat ${read_args} > seqs.fastq.gz
    pychopper -t ${task.cpus} -k ${cdna_kit} -m ${backend} ${extra} \
        seqs.fastq.gz "${meta.alias}_full_length_reads.fastq"
    awk '
        BEGIN { FS = OFS = "\\t" }
        NR == 1 {
            for (i = 1; i <= NF; i++) {
                idx[\$i] = i
            }
            print "Classification", "Value"
            next
        }
        \$idx["Category"] == "Classification" {
            print \$idx["Name"], \$idx["Value"]
        }
    ' pychopper.tsv > pychopper_summary.tsv
    mkdir -p "${meta.alias}_pychopper_output"
    find . -maxdepth 1 -mindepth 1 \
        ! -name "seqs.fastq.gz" \
        ! -name "${meta.alias}_pychopper_output" \
        -exec mv -t "${meta.alias}_pychopper_output" {} +
    """
 }
 process makeReport {
    label "wf_common"
    publishDir "${params.out_dir}", mode: 'copy', pattern: "wf-transcriptomes-report.html"
@ -101,10 +39,8 @@ process makeReport {
        tuple val(metadata), path(stats, stageAs: "stats_*")
        path "versions/*"
        path "params.json"
        path alignment_stats, stageAs: "alignment_stats/*"
        path cohort_dir, stageAs: "cohort"
        path sample_dirs, stageAs: "samples/*"
        path pychopper_dirs, stageAs: "pychopper/*"
        path sqanti_dirs, stageAs: "sqanti/*"
        path de_files
        val wf_version
@ -112,20 +48,17 @@ process makeReport {
        path "wf-transcriptomes-report.html", emit: report
    script:
        String metadata_json = new JsonBuilder(metadata).toPrettyString().replaceAll("'", "'\\\\''")
-        def report_stats = stats instanceof java.util.Collection ? stats : (stats ? [stats] : [])
+        def report_stats = (stats instanceof java.util.Collection ? stats : (stats ? [stats] : []))
-        def report_pychopper = pychopper_dirs instanceof java.util.Collection ? pychopper_dirs : (pychopper_dirs ? [pychopper_dirs] : [])
+            .findAll { it.name != OPTIONAL_FILE.name }
        String stats_args = report_stats ? "--stats ${report_stats.join(' ')}" : ""
        String pychopper_args = report_pychopper.find { it.name != OPTIONAL_FILE.name } ? "--pychopper_dir pychopper" : ""
        String de_args = de_files.name == OPTIONAL_FILE.name ? "" : "--de_dir de_analysis"
    """
    echo '${metadata_json}' > metadata.json
    workflow-glue report wf-transcriptomes-report.html \
        --metadata metadata.json \
        ${stats_args} \
        --alignment_stats_dir alignment_stats \
        --cohort_dir cohort \
        --samples_dir samples \
        ${pychopper_args} \
        --sqanti_dir sqanti \
        ${de_args} \
        --versions versions \
@ -179,7 +112,6 @@ def coerceBooleanParam(value) {
    "version",
    "igv",
    "direct_rna",
    "cdna_preprocess",
    "de_analysis",
    "analyse_unclassified",
    "analyse_fail",
@ -211,17 +143,10 @@ workflow pipeline {
        sample_sheet
        ref_genome
        ref_annotation
        pychopper_dirs
    main:
        software_versions = getVersions()
        workflow_params = getParams()
        ingress_results = collectIngressResultsInDir(
            reads.map { meta, sample_reads, stats ->
                [meta, sample_reads, stats ?: OPTIONAL_FILE]
            }
        )
        transcriptome = transcriptome_analysis(reads, ref_genome, ref_annotation, sample_sheet)
        if (params.de_analysis) {
@ -236,50 +161,53 @@ workflow pipeline {
        }
        report_input = reads
-            .map { meta, sample_reads, stats ->
+            .collect(flat: false)
-                ["all_samples", meta + [has_stats: stats as boolean], stats]
+            .map { rows ->
                def metadata = rows.collect { meta, xam, xai, stats ->
                    meta + [has_stats: stats as boolean]
                }
                def stats_rows = rows.findAll { meta, xam, xai, stats -> stats != null }
                [
                    metadata,
                    stats_rows ? stats_rows.collect { meta, xam, xai, stats -> stats } : [OPTIONAL_FILE]
                ]
            }
            .groupTuple()
            .map { report_group, metas, stats ->
                [metas, stats.findAll { it != null }]
            }
        alignment_stats = transcriptome.alignments
            .map { meta, bam, bai, flagstat -> flagstat }
            .collect()
        sample_dirs_for_report = transcriptome.sample_dirs
            .map { meta, sample_dir -> sample_dir }
            .collect()
        pychopper_dirs_for_report = pychopper_dirs
            .map { meta, pychopper_dir -> pychopper_dir }
            .ifEmpty(OPTIONAL_FILE)
            .collect()
        sqanti_dirs_for_report = transcriptome.joint_sqanti_dir
            .concat(transcriptome.sample_sqanti_dirs.map { meta, sqanti_dir -> sqanti_dir })
            .ifEmpty(OPTIONAL_FILE)
            .collect()
        // meta.src_xam is non-null if BAMs are "passed through"
        generated_alignment_outputs = reads
            .filter { meta, bam, bai, stats -> meta.src_xam == null }
            .flatMap { meta, bam, bai, stats ->
                def outdir = "cohort/alignments/${meta.alias}"
                [
                    [bam, outdir],
                    [bai, outdir],
                    [stats.resolve("bamstats.flagstat.tsv"), outdir],
                ]
            }
        report = makeReport(
            report_input,
            software_versions,
            workflow_params,
            alignment_stats,
            transcriptome.joint_dir,
            sample_dirs_for_report,
            pychopper_dirs_for_report,
            sqanti_dirs_for_report,
            de_dir,
            workflow.manifest.version
        )
        results = Channel.empty()
            .concat(ingress_results.out.map { [it, "ingress_results"] })
            .concat(report.report.map { [it, null] })
            .concat(workflow_params.map { [it, null] })
            .concat(pychopper_dirs.map { meta, pychopper_dir -> [pychopper_dir, "ingress_results/${meta.alias}"] })
            .concat(transcriptome.annotation_reference_summary.map { [it, "cohort/reference"] })
            .concat(transcriptome.unstranded_annotation.map { [it, "cohort/reference"] })
            .concat(transcriptome.joint_gtf.map { [it, "cohort"] })
@ -296,27 +224,28 @@ workflow pipeline {
            .concat(transcriptome.sample_transcript_rds.map { meta, rds -> [rds, "samples/${meta.alias}"] })
            .concat(transcriptome.sample_gene_rds.map { meta, rds -> [rds, "samples/${meta.alias}"] })
            .concat(transcriptome.sample_metadata.map { meta, metadata -> [metadata, "samples/${meta.alias}"] })
-            .concat(transcriptome.alignments.map { meta, bam, bai, flagstat -> [bam, "cohort/alignments"] })
+            .concat(generated_alignment_outputs)
            .concat(transcriptome.alignments.map { meta, bam, bai, flagstat -> [bai, "cohort/alignments"] })
            .concat(transcriptome.alignments.map { meta, bam, bai, flagstat -> [flagstat, "cohort/alignments"] })
            .concat(transcriptome.joint_sqanti_dir.map { [it, "cohort"] })
            .concat(transcriptome.sample_sqanti_dirs.map { meta, sqanti_dir -> [sqanti_dir, "samples/${meta.alias}"] })
        reference_basename = file(params.ref_genome).getName()
        if (params.igv) {
-            results = results
+            //todo ???
-                .concat(transcriptome.reference.map { [it, "igv_reference"] })
+            //results = results
-                .concat(transcriptome.reference_fai.map { [it, "igv_reference"] })
+            //    .concat(transcriptome.reference.map { [it, "igv_reference"] })
-                .concat(transcriptome.reference_gzi.map { [it, "igv_reference"] })
+            //    .concat(transcriptome.reference_fai.map { [it, "igv_reference"] })
            //    .concat(transcriptome.reference_gzi.map { [it, "igv_reference"] })
            results = Channel.empty()
-            igv_index_paths = transcriptome.reference_fai
+            //igv_index_paths = transcriptome.reference_fai
-                .map { "igv_reference/${it.getName()}" }
+            //    .map { "igv_reference/${it.getName()}" }
-                .concat(transcriptome.reference_gzi.map { "igv_reference/${it.getName()}" })
+            //    .concat(transcriptome.reference_gzi.map { "igv_reference/${it.getName()}" })
            igv_index_paths = Channel.empty()
-            igv_alignment_paths = transcriptome.alignments
+            igv_alignment_paths = reads
-                .map { meta, bam, bai, flagstat -> [
+                .map { meta, bam, bai, stat -> [
-                    "cohort/alignments/${bam.getName()}",
+                    meta.src_xam ?: "cohort/alignments/${meta.alias}/reads.bam",
-                    "cohort/alignments/${bai.getName()}"
+                    meta.src_xai ?: "cohort/alignments/${meta.alias}/reads.bam.bai"
                ] }
                .flatten()
@ -358,7 +287,7 @@ workflow {
        throw new Exception("Provide exactly one of --fastq or --bam.")
    }
    if (!params.ref_genome) {
-        throw new Exception("Provide --ref_genome.")
+        throw new Exception("Provide --ref_genome.") //todo isnt this enforced in the schema?
    }
    if (!params.ref_annotation) {
        throw new Exception("Provide --ref_annotation.")
@ -366,20 +295,20 @@ workflow {
    if (!(params.transcriptome_mode in ["discover", "fixed_annotation"])) {
        throw new Exception("--transcriptome_mode must be one of: discover, fixed_annotation.")
    }
    if (params.direct_rna && params.cdna_preprocess) {
        throw new Exception("--cdna_preprocess cannot be used together with --direct_rna.")
    }
    if (params.de_analysis && !params.sample_sheet) {
        throw new Exception("Provide --sample_sheet when running with --de_analysis.")
    }
    sample_sheet = params.sample_sheet ? file(params.sample_sheet, type: "file") : OPTIONAL_FILE
    ref_genome = file(params.ref_genome, type: "file")
    ref_annotation = file(params.ref_annotation, type: "file")
-    if (!ref_genome.exists()) {
+    prepared_reference = prepare_reference(
-        throw new Exception("--ref_genome does not exist.")
+        params.ref_genome, [
-    }
+        "output_cache": false,
        "output_mmi": false,
    ])
    ref_genome = prepared_reference.ref_tuple
    if (!ref_annotation.exists()) {
        throw new Exception("--ref_annotation does not exist.")
    }
@ -387,51 +316,31 @@ workflow {
        throw new Exception("--sample_sheet does not exist.")
    }
-    def samples
+    def ingress_args = [
        "minimap2_memory": ["31GB", "62GB"],
        "minimap2_opts": params.direct_rna ? "-ax splice -uf -k14" : "-ax splice -uf",
        "alignment_threads": 12,
        "output_xam_fmt": "bam",
        "sample": params.sample,
        "sample_sheet": params.sample_sheet,
        "analyse_unclassified": params.analyse_unclassified,
        "analyse_fail": params.analyse_fail,
        "fastcat_extra_args": "",
        "required_sample_types": [],
    ]
    if (params.fastq) {
        samples = fastq_ingress([
            "input": params.fastq,
-            "sample": params.sample,
+        ] + ingress_args, ref_genome)
            "sample_sheet": params.sample_sheet,
            "analyse_unclassified": params.analyse_unclassified,
            "analyse_fail": params.analyse_fail,
            "fastcat_extra_args": "",
            "required_sample_types": [],
            "fastq_chunk": params.fastq_chunk,
            "per_read_stats": params.wf.per_read_stats,
            "allow_multiple_basecall_models": params.wf.allow_multiple_basecall_models,
        ])
    } else {
        samples = xam_ingress([
            "input": params.bam,
-            "sample": params.sample,
+        ] + ingress_args, ref_genome)
            "sample_sheet": params.sample_sheet,
            "analyse_unclassified": params.analyse_unclassified,
            "analyse_fail": params.analyse_fail,
            "keep_unaligned": params.wf.keep_unaligned,
            "return_fastq": params.wf.return_fastq,
            "fastq_chunk": params.fastq_chunk,
            "per_read_stats": params.wf.per_read_stats,
            "allow_multiple_basecall_models": params.wf.allow_multiple_basecall_models,
        ])
    }
-    decorated_samples = samples
+    analysis_samples = samples
-        .map { meta, fname, stats -> [meta["group_key"], meta, fname, stats] }
+        .filter { meta, xam, xai, stats ->
        .groupTuple()
        .map { key, metas, fnames, statss ->
            if (fnames[0] == null) {
                fnames = null
            }
            [
                metas[0] + ["group_index": metas.collect { it["group_index"] }],
                fnames,
                statss[0]
            ]
        }
    analysis_samples = decorated_samples
        .filter { meta, sample_reads, stats ->
            if (meta.n_seqs == 0) {
                log.warn("Sample ${meta.alias} has no reads - excluded from transcriptome analysis.")
                return false
@ -442,39 +351,10 @@ workflow {
            throw new Exception("No samples with reads were available for transcriptome analysis.")
        }
-    pychopper_results = Channel.empty()
+
    processed_samples = analysis_samples
-    if (params.cdna_preprocess) {
+    pipeline_run = pipeline(processed_samples, sample_sheet, ref_genome, ref_annotation)
        grouped_samples = analysis_samples.branch { meta, sample_reads, stats ->
            to_process: sample_reads != null
            passthrough: sample_reads == null
        }
        preprocessed = preprocess_reads(
            grouped_samples.to_process
                .map { meta, sample_reads, stats -> [meta, sample_reads] }
        )
        processed_samples = grouped_samples.passthrough
            .mix(
                grouped_samples.to_process
                    .map { meta, sample_reads, stats -> [meta.alias, meta, stats] }
                    .join(preprocessed.full_len_reads)
                    .map { alias, meta, stats, full_length_reads ->
                        [meta, full_length_reads, stats]
                    }
            )
        pychopper_results = grouped_samples.to_process
            .map { meta, sample_reads, stats -> [meta.alias, meta] }
            .join(preprocessed.dir)
            .map { alias, meta, pychopper_dir ->
                [meta, pychopper_dir]
            }
    }
    pipeline_run = pipeline(processed_samples, sample_sheet, ref_genome, ref_annotation, pychopper_results)
    publishResults(pipeline_run.results)
 }
--- a/nextflow.config
+++ b/nextflow.config
@ -13,8 +13,6 @@ params {
    ref_annotation = null
    transcriptome_mode = "discover"
    direct_rna = false
    cdna_preprocess = false
    cdna_kit = "SQK-PCS109"
    de_analysis = false
    condition_column = "condition"
@ -32,8 +30,6 @@ params {
    skip_sqanti = false
    sqanti_skip_orf = true
    sqanti_extra_args = ""
    pychopper_backend = "edlib"
    pychopper_opts = ""
    aws_image_prefix = null
    aws_queue = null
--- a/nextflow_schema.json
+++ b/nextflow_schema.json
@ -96,12 +96,6 @@
          "type": "boolean",
          "default": false,
          "description": "Set this for direct RNA sequencing libraries."
        },
        "cdna_preprocess": {
          "type": "boolean",
          "default": false,
          "description": "Apply the optional cDNA preprocessing stage before alignment.",
          "help_text": "This is only relevant for cDNA libraries and must not be combined with `--direct_rna`."
        }
      },
      "required": [
@ -199,26 +193,6 @@
          "type": "number",
          "description": "Optional bambu novel discovery rate override."
        },
        "cdna_kit": {
          "type": "string",
          "default": "SQK-PCS109",
          "description": "ONT cDNA kit identifier used for pychopper preprocessing.",
          "help_text": "Only used when `--cdna_preprocess` is enabled. The workflow derives the pychopper kit code from this value."
        },
        "pychopper_backend": {
          "type": "string",
          "default": "edlib",
          "enum": [
            "edlib",
            "phmm"
          ],
          "description": "Primer-detection backend to use for pychopper preprocessing."
        },
        "pychopper_opts": {
          "type": "string",
          "default": "",
          "description": "Extra command-line options to pass to pychopper."
        },
        "skip_sqanti": {
          "type": "boolean",
          "default": false,
--- a/output_definition.json
+++ b/output_definition.json
@ -41,7 +41,7 @@
      "type": "per-sample"
    },
    "sample-bam": {
-      "filepath": "cohort/alignments/{{ alias }}.aligned.sorted.bam",
+      "filepath": "cohort/alignments/{{ alias }}/reads.bam",
      "title": "Aligned BAM",
      "description": "Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV.",
      "mime-type": "application/gzip",
@ -49,7 +49,7 @@
      "type": "per-sample"
    },
    "sample-bai": {
-      "filepath": "cohort/alignments/{{ alias }}.aligned.sorted.bam.bai",
+      "filepath": "cohort/alignments/{{ alias }}/reads.bam.bai",
      "title": "Aligned BAM index",
      "description": "Index for the aligned BAM.",
      "mime-type": "application/octet-stream",
@ -57,10 +57,10 @@
      "type": "per-sample"
    },
    "sample-flagstat": {
-      "filepath": "cohort/alignments/{{ alias }}.flagstat.txt",
+      "filepath": "cohort/alignments/{{ alias }}/bamstats.flagstat.tsv",
      "title": "Alignment summary",
-      "description": "samtools flagstat output for the aligned BAM.",
+      "description": "bamstats flagstat summary for the aligned BAM.",
-      "mime-type": "text/plain",
+      "mime-type": "text/tab-separated-values",
      "optional": false,
      "type": "per-sample"
    },
--- a/subworkflows/transcriptome.nf
+++ b/subworkflows/transcriptome.nf
@ -9,7 +9,7 @@ process prepareAnnotationReference {
    memory "6 GB"
    input:
        path ref_annotation
-        path ref_genome
+        tuple path(ref), path(ref_idx)
    output:
        stdout emit: warnings
        path "annotation.gtf", emit: annotation
@ -20,62 +20,19 @@ process prepareAnnotationReference {
    """
    workflow-glue prepare_annotation_reference \
        --annotation "${ref_annotation}" \
-        --reference "${ref_genome}" \
+        --reference "${ref}" \
        --out_dir prepared
    mv prepared/* .
    """
 }
 process buildMinimapIndex {
    label "wf_transcriptomes"
    cpus 4
    memory "32 GB"
    input:
        path reference
    output:
        tuple path("genome.mmi"), path(reference), emit: index
    script:
        String index_opts = params.minimap2_index_opts ?: ""
    """
    minimap2 -t ${task.cpus} ${index_opts} -d genome.mmi "${reference}"
    """
 }
 process alignReads {
    label "wf_transcriptomes"
    cpus { params.threads ?: 4 }
    memory "16 GB"
    input:
        tuple val(meta), path(reads), path(index), path(reference)
    output:
        tuple val(meta),
            path("${meta.alias}.aligned.sorted.bam"),
            path("${meta.alias}.aligned.sorted.bam.bai"),
            path("${meta.alias}.flagstat.txt"),
            emit: bam
    script:
        String preset = params.direct_rna ? "-ax splice -uf -k14" : "-ax splice -uf"
        String extra = params.minimap2_opts ?: ""
        def read_args = reads instanceof java.util.ArrayList ? reads.join(" ") : reads
        int sort_threads = Math.max(task.cpus as int - 1, 1)
    """
    minimap2 -t ${task.cpus} ${preset} ${extra} "${index}" ${read_args} \
        | samtools sort --write-index -@ ${sort_threads} \
            -o "${meta.alias}.aligned.sorted.bam##idx##${meta.alias}.aligned.sorted.bam.bai" -
    samtools flagstat "${meta.alias}.aligned.sorted.bam" > "${meta.alias}.flagstat.txt"
    """
 }
 process runJointBambu {
    label "wf_transcriptomes"
    cpus { params.threads ?: 4 }
    memory "32 GB"
    input:
-        path bam_files, stageAs: "bams/*"
+        tuple val(aliases), path(bams, stageAs: "bams/??.bam"), path(bais, stageAs: "bams/??.bam.bai")
        path bam_indexes, stageAs: "bams/*"
        path sample_sheet
        path annotation, stageAs: "annotation/*"
        path reference, stageAs: "reference/*"
@ -88,11 +45,16 @@ process runJointBambu {
        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 \
-        --bam_dir bams \
+        ${bams_arg} \
        ${aliases_arg} \
        ${sample_sheet_arg} \
        --annotation "${annotation}" \
        --genome "${reference}" \
@ -109,7 +71,7 @@ process runPerSampleBambu {
    cpus { params.threads ?: 4 }
    memory "24 GB"
    input:
-        tuple val(meta), path(bam), path(bai), path(flagstat)
+        tuple val(meta), path(bam), path(bai), path(stats)
        path annotation, stageAs: "annotation/*"
        path reference, stageAs: "reference/*"
    output:
@ -121,11 +83,13 @@ process runPerSampleBambu {
        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 \
-        --bam_path "${bam}" \
+        ${bams_arg} \
-        --sample_alias "${meta.alias}" \
+        ${aliases_arg} \
        --annotation "${annotation}" \
        --genome "${reference}" \
        --transcriptome_mode "${params.transcriptome_mode}" \
@ -234,40 +198,9 @@ process runPerSampleSqanti {
 }
 process faidx {
    label "wf_transcriptomes"
    cpus 1
    memory "4 GB"
    input:
        path ref
    output:
        path("${ref}.fai"), emit: fai
    script:
    """
    samtools faidx "${ref}"
    """
 }
 process gzFaidx {
    label "wf_transcriptomes"
    cpus 1
    memory "4 GB"
    errorStrategy "ignore"
    input:
        path ref
    output:
        tuple path("${ref}.fai"), path("${ref}.gzi"), emit: indexes
    script:
    """
    samtools faidx "${ref}"
    """
 }
 workflow transcriptome_analysis {
    take:
-        reads
+        alignments
        ref_genome
        ref_annotation
        sample_sheet
@ -281,22 +214,23 @@ workflow transcriptome_analysis {
        analysis_annotation = prepared_reference_annotation.annotation
        analysis_reference = prepared_reference_annotation.reference
        genome_index = buildMinimapIndex(analysis_reference)
        aligned = alignReads(
            reads.map { meta, sample_reads, stats -> [meta, sample_reads] }
                .combine(genome_index.index)
        )
        joint_bambu = runJointBambu(
-            aligned.bam.map { meta, bam, bai, flagstat -> bam }.collect(),
+            alignments
-            aligned.bam.map { meta, bam, bai, flagstat -> bai }.collect(),
+            | 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,
            analysis_annotation,
            analysis_reference
        )
-        sample_bambu = runPerSampleBambu(aligned.bam, analysis_annotation, analysis_reference)
+        sample_bambu = runPerSampleBambu(alignments, analysis_annotation, analysis_reference)
        joint_fasta = buildCohortTranscriptomeFasta(joint_bambu.gtf, analysis_reference)
        sample_fastas = buildSampleTranscriptomeFasta(sample_bambu.gtf, analysis_reference)
@ -311,22 +245,10 @@ workflow transcriptome_analysis {
            sample_sqanti_dirs = sample_sqanti.dir
        }
        compressed_reference = params.ref_genome.toLowerCase().endsWith("gz")
        if (compressed_reference) {
            ref_indexes = gzFaidx(ref_genome)
            reference_fai = ref_indexes.indexes.map { fai, gzi -> fai }
            reference_gzi = ref_indexes.indexes.map { fai, gzi -> gzi }
        } else {
            ref_indexes = faidx(ref_genome)
            reference_fai = ref_indexes.fai
            reference_gzi = Channel.empty()
        }
    emit:
        reference = ref_genome
        annotation = ref_annotation
        annotation_reference_summary = prepared_reference_annotation.summary
        unstranded_annotation = prepared_reference_annotation.unstranded
        alignments = aligned.bam
        joint_dir = joint_bambu.dir
        joint_gtf = joint_bambu.gtf
        joint_fasta = joint_fasta.fasta
@ -345,6 +267,4 @@ workflow transcriptome_analysis {
        sample_metadata = sample_bambu.transcript_metadata
        joint_sqanti_dir = joint_sqanti_dir
        sample_sqanti_dirs = sample_sqanti_dirs
        reference_fai = reference_fai
        reference_gzi = reference_gzi
 }