From 175e6c94377cc3a6e9d8025eab1865a138a77d33 Mon Sep 17 00:00:00 2001 From: Sam Nicholls Date: Tue, 26 May 2026 09:52:38 +0000 Subject: [PATCH] Per-sample pileup counts [CW-7196] --- .gitlab-ci.yml | 10 +- README.md | 31 +++- bin/infer_modkit_codes | 45 +++++ bin/mod_code_label | 44 +++++ bin/workflow_glue/check_valid_modbam.py | 36 ++++ .../tests/common/test_check_valid_modbam.py | 66 +++++++ data/mod_code_labels.tsv | 5 + docs/06_pipeline_overview.md | 28 ++- docs/07_input_parameters.md | 1 + docs/08_outputs.md | 2 + main.nf | 14 +- nextflow.config | 8 + nextflow_schema.json | 5 + output_definition.json | 16 ++ subworkflows/mods.nf | 166 ++++++++++++++++++ 15 files changed, 461 insertions(+), 16 deletions(-) create mode 100755 bin/infer_modkit_codes create mode 100755 bin/mod_code_label create mode 100644 bin/workflow_glue/check_valid_modbam.py create mode 100644 bin/workflow_glue/tests/common/test_check_valid_modbam.py create mode 100644 data/mod_code_labels.tsv create mode 100644 subworkflows/mods.nf diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index a3c38bc..261d126 100644 --- a/.gitlab-ci.yml +++ b/.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: diff --git a/README.md b/README.md index 052d99c..17f620a 100644 --- a/README.md +++ b/README.md @@ -222,7 +222,20 @@ Each sample is aligned to the supplied reference genome with `samples//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//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//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//`. -### 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//` contains alignments, independent per-sample `bambu` outputs and - optional per-sample `SQANTI3` outputs + optional per-sample `SQANTI3` outputs; `samples//mods/` is populated + when modified base tags are present in the aligned BAM + `de_analysis//` 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 | diff --git a/bin/infer_modkit_codes b/bin/infer_modkit_codes new file mode 100755 index 0000000..0028bd7 --- /dev/null +++ b/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() diff --git a/bin/mod_code_label b/bin/mod_code_label new file mode 100755 index 0000000..ddf0ddf --- /dev/null +++ b/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() diff --git a/bin/workflow_glue/check_valid_modbam.py b/bin/workflow_glue/check_valid_modbam.py new file mode 100644 index 0000000..43666af --- /dev/null +++ b/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 diff --git a/bin/workflow_glue/tests/common/test_check_valid_modbam.py b/bin/workflow_glue/tests/common/test_check_valid_modbam.py new file mode 100644 index 0000000..62ad6c8 --- /dev/null +++ b/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")) diff --git a/data/mod_code_labels.tsv b/data/mod_code_labels.tsv new file mode 100644 index 0000000..65c857a --- /dev/null +++ b/data/mod_code_labels.tsv @@ -0,0 +1,5 @@ +mod_code label +A:17596 inosine +A:a m6A +C:m m5C +T:17802 pseU diff --git a/docs/06_pipeline_overview.md b/docs/06_pipeline_overview.md index 0e73039..be1b4cb 100644 --- a/docs/06_pipeline_overview.md +++ b/docs/06_pipeline_overview.md @@ -67,7 +67,20 @@ Each sample is aligned to the supplied reference genome with `samples//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//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//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//`. -### 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//` contains alignments, independent per-sample `bambu` outputs and - optional per-sample `SQANTI3` outputs + optional per-sample `SQANTI3` outputs; `samples//mods/` is populated + when modified base tags are present in the aligned BAM + `de_analysis//` contains DE and DTU results for each contrast when differential analysis is enabled + `igv_reference/` contains the published reference indexes used for IGV diff --git a/docs/07_input_parameters.md b/docs/07_input_parameters.md index e1c877c..12105b2 100644 --- a/docs/07_input_parameters.md +++ b/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 | diff --git a/docs/08_outputs.md b/docs/08_outputs.md index d7cd4a5..a0f6e1b 100644 --- a/docs/08_outputs.md +++ b/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 | diff --git a/main.nf b/main.nf index d011df2..012bc0a 100644 --- a/main.nf +++ b/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( diff --git a/nextflow.config b/nextflow.config index 8fd3c61..f6412a5 100644 --- a/nextflow.config +++ b/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}" } diff --git a/nextflow_schema.json b/nextflow_schema.json index 84f09c1..508c1f1 100644 --- a/nextflow_schema.json +++ b/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": "", diff --git a/output_definition.json b/output_definition.json index 171da6c..2eeb5c2 100644 --- a/output_definition.json +++ b/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", diff --git a/subworkflows/mods.nf b/subworkflows/mods.nf new file mode 100644 index 0000000..8e8f67b --- /dev/null +++ b/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 +}