Per-sample pileup counts [CW-7196]

.gitlab-ci.yml

@@ -69,7 +69,7 @@ docker-run:
             - MATRIX_NAME: [
                 "int_discover_dna", "int_fixed_rna", "int_de_control_vs_control",
                 "smoke_discover", "smoke_fixed", "smoke_direct_rna", "smoke_de",
-                "mouse_de_0countquant"
+                "mouse_de_0countquant", "mods_bigwig_igv"
             ]
     rules:
         # NOTE As we're overriding the rules block for the included docker-run
@@ -181,6 +181,14 @@ docker-run:
                   --de_analysis --ref_genome ${CI_PROJECT_NAME}/data/mouse_subset_test/mouse_subset.fa \
                   --ref_annotation ${CI_PROJECT_NAME}/data/mouse_subset_test/mouse_subset.gtf.gz \
                 --direct_rna --sample_sheet ${CI_PROJECT_NAME}/data/mouse_subset_test/sample_sheet.csv"
+        # MM/ML tag test
+        - if: $MATRIX_NAME == "mods_bigwig_igv"
+          variables:
+            NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/mods_rna_subset.tar.gz -O ${CI_PROJECT_NAME}/data/mods_rna_subset.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/mods_rna_subset.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_49/gencode.v49.annotation.gtf.gz -O ${CI_PROJECT_NAME}/data/gencode.v49.annotation.gtf.gz && wget -nv https://ont-open-data.s3.amazonaws.com/references/human/GRCh38/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz -O ${CI_PROJECT_NAME}/data/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz"
+            NF_WORKFLOW_OPTS: "--bam ${CI_PROJECT_NAME}/data/mods_rna_subset/barcode01/merged.sorted.bam \
+                --ref_genome ${CI_PROJECT_NAME}/data/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz \
+                --ref_annotation ${CI_PROJECT_NAME}/data/gencode.v49.annotation.gtf.gz \
+                --igv"
 
 
 singularity-run:

README.md

@@ -222,7 +222,20 @@ Each sample is aligned to the supplied reference genome with
 `samples/<alias>/alignment/` are the main alignment files used for transcriptome
 analysis, optional [`SQANTI3`](https://github.com/conesalab/SQANTI3) QC, and optional IGV viewing.
 
-### 4. Cohort transcriptome construction
+### 4. Optional modified base summarisation
+
+When aligned BAMs contain modified base tags (`MM` and `ML`), the workflow also
+runs `modkit` on each sample alignment. It first checks which modified base
+codes are present in the BAM, then runs `modkit pileup` to produce a per-sample
+bedMethyl file and one bigWig track per requested or inferred modification
+under `samples/<alias>/mods/`.
+
+If `--mod_codes` is set, those codes are passed directly to `modkit pileup`.
+If it is omitted, the workflow infers the available `primary_base:mod_code`
+pairs from the aligned BAM with `modkit modbam check-tags`. These outputs are
+also included in the optional IGV configuration when `--igv` is enabled.
+
+### 5. 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
@@ -232,7 +245,7 @@ Before writing outputs, transcript filtering removes only transcripts with zero
 total transcript counts across samples; it does not use `fullLengthCounts` for
 this quantification filter.
 
-### 5. Independent per-sample transcriptomes
+### 6. 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
@@ -240,7 +253,7 @@ sample-specific GTF, FASTA, count tables, and metadata under
 specific transcript models without changing the shared cohort transcriptome used
 for DE/DTU.
 
-### 6. Transcript sequence generation and QC
+### 7. 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
@@ -248,7 +261,7 @@ transcriptomes and produces structural QC summaries. The cohort `SQANTI3`
 results live under `cohort/sqanti/`, while per-sample `SQANTI3`
 directories are published under `samples/<alias>/sqanti/`.
 
-### 7. Optional DE and DTU analysis
+### 8. 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
@@ -256,7 +269,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>/`.
 
-### 8. What you need to provide
+### 9. What you need to provide
 
 The workflow's analysis is controlled by a user provided genome, annotation, and
 `bambu` mode.
@@ -271,14 +284,15 @@ The workflow's analysis is controlled by a user provided genome, annotation, and
 * when `--de_analysis` is enabled, the sample sheet must contain `alias`, the
   primary condition column, and any requested columns named in `--covariates`
 
-### 9. How to read the output folder
+### 10. How to read the output folder
 
 The published outputs are organised around a small number of top-level
 directories:
 
 + `cohort/` contains the primary joint `bambu` transcriptome, count tables, and optional cohort `SQANTI3` outputs
 + `samples/<alias>/` contains alignments, independent per-sample `bambu` outputs and
-  optional per-sample `SQANTI3` outputs
+  optional per-sample `SQANTI3` outputs; `samples/<alias>/mods/` is populated
+  when modified base tags are present in the aligned BAM
 + `de_analysis/<contrast>/` contains DE and DTU results for each contrast when
   differential analysis is enabled
 + `igv_reference/` contains the published reference indexes used for IGV
@@ -340,6 +354,7 @@ directories:
 | Nextflow parameter name  | Type | Description | Help | Default |
 |--------------------------|------|-------------|------|---------|
 | threads | integer | Thread count to use for the core workflow processes. |  | 4 |
+| mod_codes | string | Comma-separated modified base codes to pass to modkit pileup. | Provide values accepted by `modkit pileup --modified-bases`, for example `A:a,C:m`. If omitted, the workflow infers `primary_base:mod_code` pairs from the BAM with `modkit modbam check-tags`. |  |
 | minimap2_opts | string | Extra command-line options to pass to minimap2. |  |  |
 | ndr | number | Optional bambu novel discovery rate override. |  |  |
 | skip_sqanti | boolean | Skip SQANTI3 transcript classification and QC. |  | False |
@@ -361,6 +376,8 @@ Output files may be aggregated including information for all samples or provided
 | Aligned BAM | samples/{{ alias }}/alignment/reads.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
 | Aligned BAM index | samples/{{ alias }}/alignment/reads.bam.bai | Index for the aligned BAM. | per-sample |
 | Alignment summary | samples/{{ alias }}/alignment/bamstats.flagstat.tsv | bamstats flagstat summary for the aligned BAM. | per-sample |
+| Modified base pileup | samples/{{ alias }}/mods/{{ alias }}.mods.bedmethyl.gz | Per-sample modkit bedMethyl pileup generated from the aligned BAM when MM and ML tags are present. | per-sample |
+| Modified base bigWig | samples/{{ alias }}/mods/{{ alias }}.mods.*.bw | Per-sample modkit bigWig tracks generated from the aligned BAM, with one file per requested or inferred modification code. | 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 |

bin/infer_modkit_codes

@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+
+import argparse
+import csv
+import sys
+
+
+def parse_args():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("modified_bases_tsv")
+    return parser.parse_args()
+
+
+def infer_codes(modified_bases_tsv):
+    codes = []
+    seen = set()
+
+    with open(modified_bases_tsv, newline="") as handle:
+        reader = csv.DictReader(handle, delimiter="\t")
+        for row in reader:
+            primary_base = row.get("primary_base", "").strip()
+            mod_code = row.get("mod_code", "").strip()
+            if not primary_base or not mod_code:
+                continue
+
+            code = f"{primary_base}:{mod_code}"
+            if code not in seen:
+                seen.add(code)
+                codes.append(code)
+
+    return codes
+
+
+def main():
+    args = parse_args()
+    codes = infer_codes(args.modified_bases_tsv)
+    if not codes:
+        sys.stderr.write("Failed to infer modified base codes from modified_bases.tsv\n")
+        raise SystemExit(1)
+
+    sys.stdout.write(",".join(codes) + "\n")
+
+
+if __name__ == "__main__":
+    main()

bin/mod_code_label

@@ -0,0 +1,44 @@
+#!/usr/bin/env python3
+
+import argparse
+import csv
+import re
+import sys
+
+
+def parse_args():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("mod_code")
+    parser.add_argument("label_tsv")
+    return parser.parse_args()
+
+
+def sanitize_label(label):
+    return re.sub(r"[^A-Za-z0-9._-]+", "_", label).strip("._-") or "mod"
+
+
+def load_labels(label_tsv):
+    labels = {}
+    with open(label_tsv, newline="") as handle:
+        reader = csv.DictReader(handle, delimiter="\t")
+        for row in reader:
+            mod_code = row.get("mod_code", "").strip()
+            label = row.get("label", "").strip()
+            if mod_code and label:
+                labels[mod_code] = label
+    return labels
+
+
+def get_mod_code_label(mod_code, labels):
+    label = labels.get(mod_code, mod_code.replace(":", "-"))
+    return sanitize_label(label)
+
+
+def main():
+    args = parse_args()
+    labels = load_labels(args.label_tsv)
+    sys.stdout.write(get_mod_code_label(args.mod_code, labels) + "\n")
+
+
+if __name__ == "__main__":
+    main()

bin/workflow_glue/check_valid_modbam.py

@@ -0,0 +1,36 @@
+#!/usr/bin/env python
+"""Check whether the input is a modbam."""
+
+import os
+import sys
+
+import pysam
+
+from .util import get_named_logger, wf_parser  # noqa: ABS101
+
+
+def main(args):
+    """Run the entry point."""
+    logger = get_named_logger("check_valid_modbam")
+    logger.info(f'Checking file: {args.bam}')
+
+    # Check the first 10K reads of the input bam for ML/MM fields
+    valid_reads = 0
+    fields = ['mm', 'ml']
+    for i, alignment in enumerate(pysam.AlignmentFile(args.bam)):
+        n_tags = len([
+            tag for (tag, val) in alignment.get_tags() if tag.lower() in fields])
+        if n_tags == 2:
+            valid_reads += 1
+            break
+        if i >= 9999:
+            break
+    if valid_reads == 0:
+        sys.exit(os.EX_DATAERR)
+
+
+def argparser():
+    """Argument parser for entrypoint."""
+    parser = wf_parser("check_valid_modbam")
+    parser.add_argument("bam", help="Input bam file")
+    return parser

bin/workflow_glue/tests/common/test_check_valid_modbam.py

@@ -0,0 +1,66 @@
+"""Tests for modified-base BAM validation."""
+
+import os
+
+import pytest
+from workflow_glue import check_valid_modbam
+
+
+class FakeAlignment:
+    """Minimal alignment stub exposing pysam's get_tags API."""
+
+    def __init__(self, tags):
+        """Store the synthetic SAM tags returned by ``get_tags``."""
+        self._tags = tags
+
+    def get_tags(self):
+        """Return the synthetic tag list for this fake alignment."""
+        return self._tags
+
+
+def _args(*argv):
+    return check_valid_modbam.argparser().parse_args(list(argv))
+
+
+def test_main_accepts_bam_with_mm_and_ml_tags(monkeypatch):
+    """A read carrying both MM and ML tags should pass validation."""
+    monkeypatch.setattr(
+        check_valid_modbam.pysam,
+        "AlignmentFile",
+        lambda _: [
+            FakeAlignment([("MM", "A+a.,0;"), ("ML", [255])]),
+        ],
+    )
+
+    check_valid_modbam.main(_args("input.bam"))
+
+
+def test_main_rejects_bam_without_mod_tags(monkeypatch):
+    """A BAM with no modified-base tags should exit with EX_DATAERR."""
+    monkeypatch.setattr(
+        check_valid_modbam.pysam,
+        "AlignmentFile",
+        lambda _: [
+            FakeAlignment([]),
+            FakeAlignment([("NM", 0)]),
+        ],
+    )
+
+    with pytest.raises(SystemExit) as exc_info:
+        check_valid_modbam.main(_args("input.bam"))
+
+    assert exc_info.value.code == os.EX_DATAERR
+
+
+def test_main_accepts_when_mod_tags_appear_later_in_scan(monkeypatch):
+    """Validation should continue scanning until it finds a tagged read."""
+    monkeypatch.setattr(
+        check_valid_modbam.pysam,
+        "AlignmentFile",
+        lambda _: [
+            FakeAlignment([("NM", 0)]),
+            FakeAlignment([("mm", "C+m.,0;"), ("ml", [200])]),
+        ],
+    )
+
+    check_valid_modbam.main(_args("input.bam"))

data/mod_code_labels.tsv

@@ -0,0 +1,5 @@
+mod_code	label
+A:17596	inosine
+A:a	m6A
+C:m	m5C
+T:17802	pseU

docs/06_pipeline_overview.md

@@ -67,7 +67,20 @@ Each sample is aligned to the supplied reference genome with
 `samples/<alias>/alignment/` are the main alignment files used for transcriptome
 analysis, optional [`SQANTI3`](https://github.com/conesalab/SQANTI3) QC, and optional IGV viewing.
 
-### 4. Cohort transcriptome construction
+### 4. Optional modified base summarisation
+
+When aligned BAMs contain modified base tags (`MM` and `ML`), the workflow also
+runs `modkit` on each sample alignment. It first checks which modified base
+codes are present in the BAM, then runs `modkit pileup` to produce a per-sample
+bedMethyl file and one bigWig track per requested or inferred modification
+under `samples/<alias>/mods/`.
+
+If `--mod_codes` is set, those codes are passed directly to `modkit pileup`.
+If it is omitted, the workflow infers the available `primary_base:mod_code`
+pairs from the aligned BAM with `modkit modbam check-tags`. These outputs are
+also included in the optional IGV configuration when `--igv` is enabled.
+
+### 5. 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
@@ -77,7 +90,7 @@ Before writing outputs, transcript filtering removes only transcripts with zero
 total transcript counts across samples; it does not use `fullLengthCounts` for
 this quantification filter.
 
-### 5. Independent per-sample transcriptomes
+### 6. 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
@@ -85,7 +98,7 @@ sample-specific GTF, FASTA, count tables, and metadata under
 specific transcript models without changing the shared cohort transcriptome used
 for DE/DTU.
 
-### 6. Transcript sequence generation and QC
+### 7. 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
@@ -93,7 +106,7 @@ transcriptomes and produces structural QC summaries. The cohort `SQANTI3`
 results live under `cohort/sqanti/`, while per-sample `SQANTI3`
 directories are published under `samples/<alias>/sqanti/`.
 
-### 7. Optional DE and DTU analysis
+### 8. 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
@@ -101,7 +114,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>/`.
 
-### 8. What you need to provide
+### 9. What you need to provide
 
 The workflow's analysis is controlled by a user provided genome, annotation, and
 `bambu` mode.
@@ -116,14 +129,15 @@ The workflow's analysis is controlled by a user provided genome, annotation, and
 * when `--de_analysis` is enabled, the sample sheet must contain `alias`, the
   primary condition column, and any requested columns named in `--covariates`
 
-### 9. How to read the output folder
+### 10. How to read the output folder
 
 The published outputs are organised around a small number of top-level
 directories:
 
 + `cohort/` contains the primary joint `bambu` transcriptome, count tables, and optional cohort `SQANTI3` outputs
 + `samples/<alias>/` contains alignments, independent per-sample `bambu` outputs and
-  optional per-sample `SQANTI3` outputs
+  optional per-sample `SQANTI3` outputs; `samples/<alias>/mods/` is populated
+  when modified base tags are present in the aligned BAM
 + `de_analysis/<contrast>/` contains DE and DTU results for each contrast when
   differential analysis is enabled
 + `igv_reference/` contains the published reference indexes used for IGV

docs/07_input_parameters.md

@@ -50,6 +50,7 @@
 | Nextflow parameter name  | Type | Description | Help | Default |
 |--------------------------|------|-------------|------|---------|
 | threads | integer | Thread count to use for the core workflow processes. |  | 4 |
+| mod_codes | string | Comma-separated modified base codes to pass to modkit pileup. | Provide values accepted by `modkit pileup --modified-bases`, for example `A:a,C:m`. If omitted, the workflow infers `primary_base:mod_code` pairs from the BAM with `modkit modbam check-tags`. |  |
 | minimap2_opts | string | Extra command-line options to pass to minimap2. |  |  |
 | ndr | number | Optional bambu novel discovery rate override. |  |  |
 | skip_sqanti | boolean | Skip SQANTI3 transcript classification and QC. |  | False |

docs/08_outputs.md

@@ -6,6 +6,8 @@ Output files may be aggregated including information for all samples or provided
 | Aligned BAM | samples/{{ alias }}/alignment/reads.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
 | Aligned BAM index | samples/{{ alias }}/alignment/reads.bam.bai | Index for the aligned BAM. | per-sample |
 | Alignment summary | samples/{{ alias }}/alignment/bamstats.flagstat.tsv | bamstats flagstat summary for the aligned BAM. | per-sample |
+| Modified base pileup | samples/{{ alias }}/mods/{{ alias }}.mods.bedmethyl.gz | Per-sample modkit bedMethyl pileup generated from the aligned BAM when MM and ML tags are present. | per-sample |
+| Modified base bigWig | samples/{{ alias }}/mods/{{ alias }}.mods.*.bw | Per-sample modkit bigWig tracks generated from the aligned BAM, with one file per requested or inferred modification code. | 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 |

main.nf

@@ -8,6 +8,8 @@ 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'
+include { mod_analysis } from './subworkflows/mods'
+
 
 
 OPTIONAL_FILE = file("$projectDir/data/OPTIONAL_FILE")
@@ -88,7 +90,6 @@ process publishResults {
     """
 }
 
-
 def coerceBooleanParam(value) {
     if (value == null || value instanceof Boolean) {
         return value
@@ -150,6 +151,7 @@ workflow pipeline {
         workflow_params = getParams()
 
         transcriptome = transcriptome_analysis(reads, ref_genome, ref_annotation, sample_sheet)
+        mods = mod_analysis(reads, ref_genome)
 
         if (params.de_analysis) {
             de_results = differential_expression(
@@ -234,6 +236,7 @@ workflow pipeline {
         }
     emit:
         results = results
+        bigwigs = mods.bigwig
 }
 
 
@@ -363,9 +366,18 @@ workflow {
             ] }
             .flatten()
 
+        // convert [alias0, [bw00...bw0N]] to [alias0, bw00] ... [aliasN, bwNN]
+        // allowing for [aliasM, bwM0] if only one bw is output because ... nextflow
+        // and use the anticipated output location
+        igv_bigwigs = pipeline_run.bigwigs
+            .flatMap { alias, paths ->
+                (paths instanceof List ? paths : [paths]).collect { path -> "${alias},samples/${alias}/mods/${path.name}" }
+            }
+
         igv_files = igv_files
             .concat(igv_index_paths)
             .concat(igv_alignment_paths)
+            .concat(igv_bigwigs)
             .collectFile(name: "igv-files.txt", newLine: true, sort: false)
 
         igv_conf = configure_igv(

nextflow.config

@@ -24,6 +24,7 @@ params {
     fastq_chunk = null
 
     threads = 4
+    mod_codes = null
     minimap2_index_opts = ""
     minimap2_opts = ""
     ndr = null
@@ -57,6 +58,7 @@ params {
         ]
         common_sha = "sha21d552f9910c575766e5d465fcb7b52fefda4b79"
         container_sha = "sha02e44f706d88fa29d8344b78479f187db7eec4ec"
+        modkit_sha = "sha489d708a48c66368e5d1e118538e5dca68203a64"
         sqanti_sha = "sha5bd775836492699e2537ebf846098eb117191d87"
         agent = null
         epi2me_instance = null
@@ -81,6 +83,9 @@ process {
     withLabel:wf_transcriptomes {
         container = "ontresearch/wf-transcriptomes-core:${params.wf.container_sha}"
     }
+    withLabel:modkit {
+        container = "ontresearch/modkit:${params.wf.modkit_sha}"
+    }
     withLabel:wf_transcriptomes_sqanti {
         container = "ontresearch/wf-transcriptomes-sqanti:${params.wf.sqanti_sha}"
     }
@@ -117,6 +122,9 @@ profiles {
             withLabel:wf_transcriptomes {
                 container = "${params.aws_image_prefix}-wf-transcriptomes-core:${params.wf.container_sha}"
             }
+            withLabel:modkit {
+                container = "${params.aws_image_prefix}-modkit:${params.wf.modkit_sha}"
+            }
             withLabel:wf_transcriptomes_sqanti {
                 container = "${params.aws_image_prefix}-wf-transcriptomes-sqanti:${params.wf.sqanti_sha}"
             }

nextflow_schema.json

@@ -178,6 +178,11 @@
           "default": 4,
           "description": "Thread count to use for the core workflow processes."
         },
+        "mod_codes": {
+          "type": "string",
+          "description": "Comma-separated modified base codes to pass to modkit pileup.",
+          "help_text": "Provide values accepted by `modkit pileup --modified-bases`, for example `A:a,C:m`. If omitted, the workflow infers `primary_base:mod_code` pairs from the BAM with `modkit modbam check-tags`."
+        },
         "minimap2_index_opts": {
           "type": "string",
           "default": "",

output_definition.json

@@ -32,6 +32,22 @@
       "optional": false,
       "type": "per-sample"
     },
+    "sample-bedmethyl": {
+      "filepath": "samples/{{ alias }}/mods/{{ alias }}.mods.bedmethyl.gz",
+      "title": "Modified base pileup",
+      "description": "Per-sample modkit bedMethyl pileup generated from the aligned BAM when MM and ML tags are present.",
+      "mime-type": "application/gzip",
+      "optional": true,
+      "type": "per-sample"
+    },
+    "sample-mod-bigwig": {
+      "filepath": "samples/{{ alias }}/mods/{{ alias }}.mods.*.bw",
+      "title": "Modified base bigWig",
+      "description": "Per-sample modkit bigWig tracks generated from the aligned BAM, with one file per requested or inferred modification code.",
+      "mime-type": "application/octet-stream",
+      "optional": true,
+      "type": "per-sample"
+    },
     "annotation-reference-summary": {
       "filepath": "cohort/reference/annotation_reference_summary.json",
       "title": "Reference and annotation preparation summary",

subworkflows/mods.nf

@@ -0,0 +1,166 @@
+// Check that the bam has modifications
+process validate_modbam {
+    label "wf_common"
+    cpus 1
+    memory 4.GB
+    input:
+        tuple val(meta),
+            path(alignment),
+            path(alignment_index),
+            val(alignment_stats)
+    output:
+        tuple val(meta),
+            path(alignment),
+            path(alignment_index),
+            val(alignment_stats),
+            env(valid)
+
+    script:
+    """
+    valid=0
+    workflow-glue check_valid_modbam ${alignment} || valid=\$?
+
+    # Allow EX_OK and EX_DATAERR, otherwise explode
+    if [ \$valid -ne 0 ] && [ \$valid -ne 65 ]; then
+        exit 1
+    fi
+    """
+}
+
+process runModkitPileup {
+    label "modkit"
+    cpus { params.threads ?: 4 }
+    memory "16 GB"
+    input:
+        tuple val(alias),
+            val(meta),
+            path(alignment),
+            path(alignment_index),
+            val(alignment_stats),
+            val(mod_codes)
+        tuple path(reference),
+            path(reference_index)
+    output:
+        tuple val(alias),
+            path("${alias}.mods.bedmethyl.gz"),
+            emit: bedmethyl
+    publishDir "${params.out_dir}/${output_key}/mods"
+
+    script:
+    output_key = alias == "cohort" ? "cohort" : "samples/${alias}"  // nodef
+    String modified_bases_arg = "--modified-bases " + mod_codes
+        .split(',')
+        .join(' ')
+    """
+    modkit pileup \
+        "${alignment}" \
+        "${alias}.mods.bedmethyl.gz" \
+        ${modified_bases_arg} \
+        --reference "${reference}" \
+        --threads ${task.cpus} \
+        --bgzf
+    """
+}
+
+process modkit_tobigwig {
+    label "modkit"
+    cpus 4
+    memory "2 GB"
+    input:
+        tuple path(reference),
+            path(reference_index)
+        tuple val(alias),
+            path(bedmethyl),
+            val(mod_codes)
+        path mod_code_labels
+    output:
+        tuple val(alias),
+            path("${alias}.mods.*.bw"),
+            emit: bigwig
+    publishDir "${params.out_dir}/${output_key}/mods"
+
+    script:
+    output_key = alias == "cohort" ? "cohort" : "samples/${alias}"  // nodef
+    String mod_code_args = mod_codes
+        .split(',')
+        .join(' ')
+    """
+    for mod_code in ${mod_code_args}; do
+        mod_code_value="\${mod_code#*:}"
+        mod_label=\$(mod_code_label "\${mod_code}" "${mod_code_labels}")
+        zcat "${bedmethyl}" | \
+            modkit bedmethyl tobigwig \
+                --sizes "${reference_index}" \
+                --nthreads ${task.cpus} \
+                --mod-codes "\${mod_code_value}" \
+                - \
+                "${alias}.mods.\${mod_label}.bw"
+    done
+    """
+}
+
+process inferModkitBases {
+    label "modkit"
+    cpus 1
+    memory "4 GB"
+    input:
+        tuple val(alias),
+            val(meta),
+            path(alignment),
+            path(alignment_index),
+            val(alignment_stats)
+    output:
+        tuple val(alias),
+            env(mod_codes)
+
+    script:
+    """
+    modkit modbam check-tags "${alignment}" --num-reads 10000 --mapped-only --out-dir check_tags
+    infer_modkit_codes check_tags/modified_bases.tsv > mod_codes.txt
+    if [ ! -s mod_codes.txt ]; then
+        echo "Failed to infer modified base codes from ${alignment}" >&2
+        exit 1
+    fi
+
+    read -r mod_codes < mod_codes.txt
+    """
+}
+
+workflow mod_analysis {
+    take:
+        xams
+        ref_genome
+    main:
+        // Check inputs have modtags, we'll early abort mod analysis for any that don't
+        validate_modbam(xams)
+            .branch {
+                nomods: it[-1] == '65'
+                    return it[0].alias
+                mods: it[-1] == '0'
+                    return [it[0].alias] + it[0..-2]  // prepend alias for joining and drop exit_code marker
+            }
+            .set{xams_with}
+
+        // warn for samples without mods
+        xams_with.nomods.subscribe {
+            log.warn "Input ${it} does not contain modified base tags. Was a modified basecalling model selected when basecalling this data?"
+        }
+
+        // determine what mods to ask modkit pileup for
+        sample_modcodes = params.mod_codes
+            ? xams_with.mods.map { [it[0], params.mod_codes.trim()] }  // cross all aliases with user mod_codes
+            : inferModkitBases(xams_with.mods)                         // otherwise infer per-sample from modbam
+        mod_samples = xams_with.mods.join(sample_modcodes)
+
+        pileup = runModkitPileup(mod_samples, ref_genome)
+
+        sample_bigwigs = modkit_tobigwig(
+            ref_genome,
+            pileup.bedmethyl.join(sample_modcodes),
+            file("$projectDir/data/mod_code_labels.tsv")
+        )
+
+    emit:
+        bedmethyl = pileup.bedmethyl
+        bigwig = sample_bigwigs.bigwig
+}