From c21588761bde14da300e9ecf9e5d25c3e1122a8e Mon Sep 17 00:00:00 2001 From: Sam Nicholls Date: Tue, 9 Jun 2026 11:58:17 +0000 Subject: [PATCH] wf-transcriptomes v2.0.0 --- .gitlab-ci.yml | 27 +++- CHANGELOG.md | 6 +- README.md | 19 ++- bin/workflow_glue/check_valid_modbam.py | 36 ------ bin/workflow_glue/report.py | 59 +++++---- .../tests/common/test_check_sample_sheet.py | 44 ++++--- .../tests/common/test_check_valid_modbam.py | 66 ---------- bin/workflow_glue/volcano.py | 82 +++++++++--- .../wfg_helpers/check_bam_headers_in_dir.py | 121 +++++++++++++----- .../wfg_helpers/validators/wf.py | 6 +- bin/workflow_glue_r/R/common.R | 21 --- bin/workflow_glue_r/R/de_analysis.R | 4 - .../tests/testthat/test_de_analysis.R | 33 ----- docs/04_install_and_run.md | 18 +++ docs/07_input_parameters.md | 1 - evaluation/tests.sh | 60 --------- lib/common.nf | 2 +- lib/ingress.nf | 57 +++++++-- main.nf | 99 +++++++++++++- nextflow.config | 5 +- nextflow_schema.json | 7 +- subworkflows/differential_expression.nf | 2 +- subworkflows/mods.nf | 46 +------ subworkflows/transcriptome.nf | 4 +- 24 files changed, 440 insertions(+), 385 deletions(-) delete mode 100644 bin/workflow_glue/check_valid_modbam.py delete mode 100644 bin/workflow_glue/tests/common/test_check_valid_modbam.py delete mode 100755 evaluation/tests.sh diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d1bc4ca..f6db86e 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -16,9 +16,7 @@ variables: PYTEST_TESTS_PATH: "bin/workflow_glue/tests/common" RTEST_CONTAINER_NAME: "wf-transcriptomes-core" RTEST_CONTAINER_CONFIG_KEY: "container_sha" - WF_TEMPLATE_ENFORCEMENT_BRANCH: "ruo-next" - RELEASE_WORKFLOW: "no" # do not release this workflow! - RELEASE_WITHOUT_DEMO: "yes" + WF_TEMPLATE_ENFORCEMENT_REF: "v6.1.0-rc1" pytest_wfcontainer: extends: pytest @@ -74,7 +72,8 @@ 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", "mods_bigwig_igv", "mismatch-sample-alias" + "mouse_de_0countquant", "mods_bigwig_igv", "mismatch-sample-alias", + "mouse_splice_error", "mouse_splice_realign" ] rules: # NOTE As we're overriding the rules block for the included docker-run @@ -202,6 +201,26 @@ docker-run: --ref_annotation ${CI_PROJECT_NAME}/data/gencode.v22.annotation.chr20.gtf \ --sample_sheet ${CI_PROJECT_NAME}/data/mods_rna_subset/sample_sheet.csv \ --igv" + - if: $MATRIX_NAME == "mouse_splice_error" + variables: + NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/mouse_splice_fail.tar.gz -O ${CI_PROJECT_NAME}/data/mouse_splice_fail.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/mouse_splice_fail.tar.gz -C ${CI_PROJECT_NAME}/data/" + NF_WORKFLOW_OPTS: "--bam ${CI_PROJECT_NAME}/data/mouse_splice_fail/samples \ + --de_analysis --ref_genome ${CI_PROJECT_NAME}/data/mouse_splice_fail/mouse_subset.fa \ + --ref_annotation ${CI_PROJECT_NAME}/data/mouse_splice_fail/mouse_subset.gtf.gz \ + --sample_sheet ${CI_PROJECT_NAME}/data/mouse_splice_fail/sample_sheet.csv" + ASSERT_NEXTFLOW_FAILURE: "1" + AFTER_NEXTFLOW_CMD: > + grep -qF 'Cannot proceed with mixed splice-aware CIGAR evidence in input BAMs.' .nextflow.log + - if: $MATRIX_NAME == "mouse_splice_realign" + variables: + NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/mouse_splice_realign.tar.gz -O ${CI_PROJECT_NAME}/data/mouse_splice_realign.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/mouse_splice_realign.tar.gz -C ${CI_PROJECT_NAME}/data/" + NF_WORKFLOW_OPTS: "--bam ${CI_PROJECT_NAME}/data/mouse_splice_realign/samples \ + --de_analysis --ref_genome ${CI_PROJECT_NAME}/data/mouse_splice_realign/mouse_subset.fa \ + --ref_annotation ${CI_PROJECT_NAME}/data/mouse_splice_realign/mouse_subset.gtf.gz \ + --sample_sheet ${CI_PROJECT_NAME}/data/mouse_splice_realign/sample_sheet.csv" + AFTER_NEXTFLOW_CMD: > + grep -q 'bam_alignment' ${CI_PROJECT_NAME}/execution/trace.txt && grep -q 'No input BAMs appear to contain splice-aware CIGAR evidence. The workflow will realign all inputs.' .nextflow.log + singularity-run: diff --git a/CHANGELOG.md b/CHANGELOG.md index 1b4e2a0..06894ee 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,8 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [v2.0.0-rc1] + +## [v2.0.0] This release refreshes `wf-transcriptomes` around a new reference-guided transcriptomics workflow built on `bambu`, with `SQANTI3` transcript classification and QC, `DESeq2` for differential gene expression, `DEXSeq` for differential transcript usage, and per-sample modified base summarisation with `modkit` when modification tags are present in aligned BAMs. @@ -32,9 +33,10 @@ This release refreshes `wf-transcriptomes` around a new reference-guided transcr - Sample-level hierarchical clustering, PCA and distance heatmap plots - Contrast-level interactive volcano plots ### Removed +- Dependence on the older `StringTie/GffCompare/Salmon` transcriptomics pathway. - `--transcriptome_source` parameter; use `--transcriptome_mode` instead. - `--ref_transcriptome` parameter; use `--transcriptome_mode fixed_annotation` together with `--ref_genome` and `--ref_annotation`. -- Dependence on the older `StringTie/GffCompare/Salmon` transcriptomics pathway. +- `--threads` parameter; advanced users may use Nextflow process selectors to override the preset per-process CPU and memory limits. ## [v1.7.2] This patch release of wf-transcriptomes updates internal workflow naming, and does not affect any workflow outputs. diff --git a/README.md b/README.md index c16605a..be7d1dd 100644 --- a/README.md +++ b/README.md @@ -112,6 +112,24 @@ the following command: nextflow pull epi2me-labs/wf-transcriptomes ``` +A demo dataset is provided for testing of the workflow. +It can be downloaded and unpacked using the following commands: +``` +wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/v2/wf-transcriptomes-demo.tar.gz +tar -xzvf wf-transcriptomes-demo.tar.gz +``` +The workflow can then be run with the downloaded demo data using: +``` +nextflow run epi2me-labs/wf-transcriptomes \ + --bam 'wf-transcriptomes-demo/samples' \ + --de_analysis \ + --direct_rna \ + --ref_annotation 'wf-transcriptomes-demo/gencode.v22.annotation.chr20.gtf' \ + --ref_genome 'wf-transcriptomes-demo/hg38_chr20.fa' \ + --sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv' \ + -profile standard +``` + @@ -351,7 +369,6 @@ The workflow's analysis is controlled by a user provided genome, annotation, and | 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`. | | | force_alignment | boolean | Force re-alignment of input BAM files. | Read alignment is skipped if the existing sequence names in the aligned BAM match the provided reference. Enable this if the existing alignments used incorrect minimap2 presets (e.g. missing --splice or direct RNA settings). | False | | ndr | number | Optional bambu novel discovery rate override. | Lower values are more conservative (higher precision), while higher values are more permissive (higher novel-discovery sensitivity). See the [`bambu` repository](https://github.com/GoekeLab/bambu) for method details. | | diff --git a/bin/workflow_glue/check_valid_modbam.py b/bin/workflow_glue/check_valid_modbam.py deleted file mode 100644 index 43666af..0000000 --- a/bin/workflow_glue/check_valid_modbam.py +++ /dev/null @@ -1,36 +0,0 @@ -#!/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/report.py b/bin/workflow_glue/report.py index fabec8b..ddc7ea4 100644 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -386,21 +386,28 @@ def _mod_summary_matrix_style(): """ -def _contrast_results(de_dir, filename, n=None): +def _contrast_results(de_dir, filename): """Return a dict of per-contrast result DataFrames read from filename.""" tables = {} for contrast_dir in sorted(Path(de_dir).iterdir()): if not contrast_dir.is_dir(): continue - table = _read_table(contrast_dir / filename) + # Enforce str dtype in case of all Nan values. + table = _read_table( + contrast_dir / filename, + dtype={'gene_name': str, 'transcript_name': str} + ) if table is None or table.empty: continue - data = table - if n is not None: - data = data.head(n) - if "padj" in data.columns: - data.sort_values("padj", ascending=True, inplace=True) - tables[contrast_dir.name] = data + + if 'gene_name' in table.columns: + table["gene_name"] = table["gene_name"].fillna("-") + if 'transcript_name' in table.columns: + table["transcript_name"] = table["transcript_name"].fillna("-") + + if "padj" in table.columns: + table.sort_values("padj", ascending=True, inplace=True) + tables[contrast_dir.name] = table return tables @@ -548,17 +555,23 @@ def _heatmap_style(): def _volcano_style(): return """ - .volcano-table-grid { + .volcano-plot-grid { display: grid; - grid-template-columns: repeat(2, minmax(0, 1fr)); - gap: 20px 10px; + grid-template-columns: minmax(0, 1fr) minmax(430px, 34%); + gap: 18px; align-items: start; } - .volcano-table-grid > * { + .volcano-plot-grid > *, + .volcano-side-panel > * { min-width: 0; } - @media screen and (max-width: 1000px) { - .volcano-table-grid { + .volcano-side-panel { + display: grid; + gap: 10px; + align-items: start; + } + @media screen and (max-width: 1100px) { + .volcano-plot-grid { grid-template-columns: 1fr; } } @@ -1417,10 +1430,11 @@ def main(args): h3("Gene expression volcano Plot") gn_vol, gn_class_table, gn_selected_table = volcano(table) - EZChart(gn_vol, width="100%", height="550") - with div(_class="volcano-table-grid"): - EZChart(gn_class_table, width="100%", height="auto") - EZChart(gn_selected_table, width="100%", height="auto") + with div(_class="volcano-plot-grid"): + EZChart(gn_vol, width="100%", height="550") + with div(_class="volcano-side-panel"): + EZChart(gn_class_table, width="100%", height="auto") + EZChart(gn_selected_table, width="100%", height="auto") with report.add_section("Differential transcript usage", "DTU"): p( @@ -1505,10 +1519,11 @@ def main(args): h3("Transcript expression volcano Plot") tr_vol, tr_class_table, tr_selected_table = volcano(dtu_table) - EZChart(tr_vol, width="100%", height="550") - with div(_class="volcano-table-grid"): - EZChart(tr_class_table, width="100%", height="auto") - EZChart(tr_selected_table, width="100%", height="auto") + with div(_class="volcano-plot-grid"): + EZChart(tr_vol, width="100%", height="550") + with div(_class="volcano-side-panel"): + EZChart(tr_class_table, width="100%", height="auto") + EZChart(tr_selected_table, width="100%", height="auto") else: p("No DTU results available for this contrast.") diff --git a/bin/workflow_glue/tests/common/test_check_sample_sheet.py b/bin/workflow_glue/tests/common/test_check_sample_sheet.py index ff9721e..bba24e0 100644 --- a/bin/workflow_glue/tests/common/test_check_sample_sheet.py +++ b/bin/workflow_glue/tests/common/test_check_sample_sheet.py @@ -26,13 +26,12 @@ def _run_check_sample_sheet(sample_sheet_path, params_json, capsys): [ [ "condition", "batch_1", "site.2", - "S01", "a1", "sample_group", - "01", "1batch", "2.site-3" + "S01", "a1", "sample_group" ], ], ) def test_validate_r_formula_names_accepts_valid_names(names): - """R formula names accept alnum-led names with safe punctuation.""" + """R formula names accept letter-led names with safe punctuation.""" _validate_r_formula_names(names) @@ -45,6 +44,9 @@ def test_validate_r_formula_names_accepts_valid_names(names): ["size:batch"], ["size|batch"], ["sueño"], + ["01"], + ["1batch"], + ["2.site-3"], [""] ], ) @@ -240,7 +242,7 @@ def test_check_sample_sheet_rejects_unsafe_covariate_values( """Configured covariate values must use safe contrast/design characters.""" sample_sheet_path = tmp_path / "sample_sheet_invalid_covariate_value.csv" sample_sheet_path.write_text( - "sample_name,condition,good-batch,site\n" + "sample_name,condition,good_batch,site\n" "test_name1,control,b1,s1\n" "test_name2,control,b2,s2\n" "test_name3,treated,b 1,s1\n" @@ -252,7 +254,7 @@ def test_check_sample_sheet_rejects_unsafe_covariate_values( json.dumps({ "de_analysis": True, "condition_column": "condition", - "covariates": "good-batch,site", + "covariates": "good_batch,site", }) ) @@ -260,16 +262,16 @@ def test_check_sample_sheet_rejects_unsafe_covariate_values( assert out.startswith( "Covariate value names must be safe for R formulas. " - "Invalid names: b 1. Names must start with a letter or number and " - "contain only letters, numbers, underscores, dots, and hyphens. " - "(column: good-batch, line: 3)" + "Invalid names: b 1. Names must start with a letter and " + "contain only letters, numbers, underscores, and dots. " + "(column: good_batch, line: 3)" ) def test_check_sample_sheet_rejects_unsafe_column_names( tmp_path, capsys ): - """Configured DE condition and covariate column names may include hyphens.""" + """Configured DE design names may not include hyphens or spaces.""" sample_sheet_path = tmp_path / "sample_sheet_invalid_design_name.csv" sample_sheet_path.write_text( "sample_name,condition-column,bad covariate\n" @@ -289,12 +291,14 @@ def test_check_sample_sheet_rejects_unsafe_column_names( ) out = _run_check_sample_sheet(sample_sheet_path, params_json, capsys) - assert out.startswith( + assert out.splitlines() == [ + "Design column names must be safe for R formulas. " + "Invalid names: condition-column. Names must start with a letter and " + "contain only letters, numbers, underscores, and dots.", "Covariate column names must be safe for R formulas. " - "Invalid names: bad covariate. Names must start with a letter or " - "number and contain only letters, numbers, underscores, dots, and " - "hyphens." - ) + "Invalid names: bad covariate. Names must start with a letter and " + "contain only letters, numbers, underscores, and dots.", + ] def test_check_sample_sheet_rejects_unsafe_condition_values( @@ -323,13 +327,13 @@ def test_check_sample_sheet_rejects_unsafe_condition_values( assert out.splitlines() == [ "Condition value names must be safe for R formulas. " - "Invalid names: control/group. Names must start with a letter or " - "number and contain only letters, numbers, underscores, dots, and " - "hyphens. (column: condition, line: 1)", + "Invalid names: control/group. Names must start with a letter and " + "contain only letters, numbers, underscores, and dots. " + "(column: condition, line: 1)", "Condition value names must be safe for R formulas. " - "Invalid names: control/group. Names must start with a letter or " - "number and contain only letters, numbers, underscores, dots, and " - "hyphens. (column: condition, line: 2)", + "Invalid names: control/group. Names must start with a letter and " + "contain only letters, numbers, underscores, and dots. " + "(column: condition, line: 2)", ] diff --git a/bin/workflow_glue/tests/common/test_check_valid_modbam.py b/bin/workflow_glue/tests/common/test_check_valid_modbam.py deleted file mode 100644 index 62ad6c8..0000000 --- a/bin/workflow_glue/tests/common/test_check_valid_modbam.py +++ /dev/null @@ -1,66 +0,0 @@ -"""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/bin/workflow_glue/volcano.py b/bin/workflow_glue/volcano.py index 7d70945..e2e0a6d 100644 --- a/bin/workflow_glue/volcano.py +++ b/bin/workflow_glue/volcano.py @@ -1,6 +1,6 @@ """Interactive volcano plot visualization for differential expression analysis.""" -from bokeh.events import Tap +from bokeh.events import DocumentReady, Tap from bokeh.layouts import column as bokeh_column, row as bokeh_row from bokeh.models import ( AutocompleteInput, @@ -102,15 +102,7 @@ def _volcano_source_data(data, fold_threshold=1, p_threshold=0.05): if not has_gene_name: data["gene_name"] = "" - data["gene_group"] = data["GENEID"] if "GENEID" in data.columns else "" - if is_transcript_plot: - if "gene_name" in data.columns: - named_group = data["gene_name"].astype(bool) - data.loc[named_group, "gene_group"] = data.loc[named_group, "gene_name"] - else: - named_group = data["GENEID"].astype(bool) - data.loc[named_group, "gene_group"] = data.loc[named_group, "GENEID"] - + data["gene_group"] = data["GENEID"] data["selected_label"] = "" # Transformations and thresholds @@ -306,6 +298,25 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): _pos_means = source_data["mean_expression"][source_data["mean_expression"] > 0] mean_min_log = float(_pos_means.min()) if len(_pos_means) > 0 else 0.01 y_threshold = -np.log10(p_threshold) + responsive_slider_title_stylesheet = """ + @media screen and (max-width: 1100px) { + .bk-slider-title { + display: block; + max-width: 100%; + overflow: hidden; + text-overflow: ellipsis; + white-space: nowrap; + } + .bk-slider-title:hover { + overflow: visible; + text-overflow: clip; + white-space: normal; + position: relative; + z-index: 1; + background: white; + } + } + """ source = ColumnDataSource(source_data) selection_state = ColumnDataSource({"indices": [[]]}) @@ -577,8 +588,9 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): source=selected_source, columns=table_columns, autosize_mode="force_fit", - height=310, + height=355, index_position=None, + visible=False, ) legend_table = DataTable( source=legend_source, @@ -654,6 +666,7 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): step=0.1, show_value=False, sizing_mode="stretch_width", + stylesheets=[responsive_slider_title_stylesheet], ) fold_input = TextInput( title="", @@ -668,6 +681,7 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): step=0.1, show_value=False, sizing_mode="stretch_width", + stylesheets=[responsive_slider_title_stylesheet], ) p_input = TextInput( title="", @@ -681,6 +695,7 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): selection_state=selection_state, highlight_source=highlight_source, selected_source=selected_source, + selected_table=selected_table, legend_source=legend_source, fold_slider=fold_slider, fold_input=fold_input, @@ -783,6 +798,8 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): selectedData.volcano_class.push(data.volcano_class[index]); }); selected_source.data = selectedData; + selected_table.visible = selected.length > 0; + selected_table.height = window.innerWidth <= 1100 ? 180 : 355; if (view_toggle.active) { const highlightData = { log2FoldChange: [], @@ -1060,7 +1077,8 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): "Search by gene/transcript…" if is_transcript_plot else "Search by gene…"), completions=[], min_characters=1, - width=240, + height=32, + width=225, ) search_input.js_on_change("value_input", CustomJS( args=dict( @@ -1169,7 +1187,8 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): clear_button = Button( label="Clear selection", button_type="default", - width=130, + height=32, + width=120, ) clear_button.js_on_click(CustomJS( args=dict( @@ -1186,6 +1205,7 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): export_button = Button( label="\u2B07 TSV", button_type="light", + height=32, width=55 ) @@ -1232,8 +1252,15 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): controls = bokeh_row( fold_slider, fold_input, p_slider, p_input, sizing_mode="stretch_width") toggle_row = bokeh_row( - plot_mode_toggle, view_toggle, - gene_select_toggle, select_significant_button + plot_mode_toggle, + view_toggle, + gene_select_toggle, + select_significant_button, + sizing_mode="stretch_width", + styles={ + "flex-wrap": "wrap", + "row-gap": "6px", + }, ) volcano_ma_plot._fig = bokeh_column( bokeh_row(controls, sizing_mode="stretch_width"), @@ -1255,8 +1282,6 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): sizing_mode="stretch_width", styles={ "padding-top": "4px", - "padding-left": "35px", - "padding-right": "20px" }, ) selected_plot = BokehPlot() @@ -1265,8 +1290,27 @@ def volcano(data, fold_threshold=1, p_threshold=0.05): bokeh_row( search_input, clear_button, export_button) ), - selected_table, sizing_mode="stretch_width", + selected_table, + sizing_mode="stretch_width", styles={ - "padding-top": "25px", + "padding-top": "8px", }) + selected_plot._fig.js_on_event(DocumentReady, CustomJS( + args=dict( + selected_source=selected_source, + selected_table=selected_table, + ), + code=""" + const listenerKey = "volcano_selected_table_resize_" + selected_source.id; + function updateSelectedTableLayout() { + selected_table.visible = selected_source.data.source_index.length > 0; + selected_table.height = window.innerWidth <= 1100 ? 180 : 355; + } + updateSelectedTableLayout(); + if (!window[listenerKey]) { + window[listenerKey] = true; + window.addEventListener("resize", updateSelectedTableLayout); + } + """, + )) return volcano_ma_plot, classes_table, selected_plot diff --git a/bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py b/bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py index b1cbabc..2b44ca1 100644 --- a/bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py +++ b/bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py @@ -8,37 +8,69 @@ import pysam from ..util import get_named_logger, wf_parser # noqa: ABS101 -def extract_header_info(xam_file, has_reads, check_ref=False): +READ_SCAN_LIMIT = 10000 + + +def alignment_has_cigar_op(alignment, cigar_op): + """Return whether an alignment has a splice/reference-skip CIGAR op.""" + return any(op == cigar_op for op, _ in (alignment.cigartuples or [])) + + +def alignment_has_tags(alignment, tags): + """Return whether an alignment carries all tags in a list of SAM tags.""" + alignment_tags = {tag.upper() for tag, _ in alignment.get_tags()} + return all(tag in alignment_tags for tag in tags) + + +def check_n_reads( + alignment_file, + check_modbase_tags=False, + check_splice_cigars=False, + read_scan_limit=READ_SCAN_LIMIT, +): + """Extract bounded read-level facts from an open XAM file.""" + read_facts = { + "has_reads": False, + } + if check_splice_cigars: + read_facts["has_splice_cigars"] = False + if check_modbase_tags: + read_facts["has_modbase_tags"] = False + + for i, alignment in enumerate(alignment_file.fetch(until_eof=True)): + if i >= read_scan_limit: + break + read_facts["has_reads"] = True + if check_splice_cigars and not read_facts["has_splice_cigars"]: + read_facts["has_splice_cigars"] = alignment_has_cigar_op( + alignment, cigar_op=3) + if check_modbase_tags and not read_facts["has_modbase_tags"]: + read_facts["has_modbase_tags"] = alignment_has_tags( + alignment, tags=["MM", "ML"]) + + if all(read_facts.values()): + break + + return read_facts + + +def check_header(alignment_file, 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 + # 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 alignment_file.header.get("SQ", [])] - 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 = alignment_file.header.get("HD") - hd_lines = f.header.get("HD") + xam_reflen = None + if check_ref: + xam_reflen = set(zip(alignment_file.references, alignment_file.lengths)) - 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 + return sq_lines, hd_lines, xam_reflen def compare_ref_lengths(xam_reflen, ref_reflen, logger): @@ -107,16 +139,37 @@ def main(args): requires_realign = False any_sq_lines = False # Track if any file has @SQ lines has_reads = False + # a bit hacky to check these specific XAM properties here but this is a very good + # place to do it. in future i'd like to abstract this out to be more like how sample + # sheets are checked, with delegation to the workflow scope for determining checks + has_splice_cigars = False + has_modbase_tags = False for xam_file in target_files: - sq_lines, hd_lines, xam_reflen, has_reads = extract_header_info( - xam_file, has_reads, check_ref=bool(args.ref)) - - if sq_lines is None: + try: + xam_fh = pysam.AlignmentFile(xam_file, check_sq=False) + except (ValueError, IOError): # File couldn't be opened logger.error(f"Failed to open {xam_file}") continue + with xam_fh: + sq_lines, hd_lines, xam_reflen = check_header( + xam_fh, + check_ref=bool(args.ref), + ) + if not (has_reads and has_splice_cigars and has_modbase_tags): + read_facts = check_n_reads( + xam_fh, + check_splice_cigars=not has_splice_cigars, + check_modbase_tags=not has_modbase_tags, + ) + has_reads = has_reads or read_facts["has_reads"] + has_splice_cigars = ( + has_splice_cigars or read_facts["has_splice_cigars"]) + has_modbase_tags = ( + has_modbase_tags or read_facts["has_modbase_tags"]) + if sq_lines: any_sq_lines = True @@ -157,9 +210,11 @@ def main(args): f"IS_UNALIGNED={int(is_unaligned)};" + f"MIXED_SQ_HEADERS={int(mixed_sq_headers)};" + f"IS_SORTED={int(sorted_xam)};" + - f"HAS_READS={int(has_reads)};" + f"HAS_READS={int(has_reads)};" + + f"HAS_SPLICE_CIGARS={int(has_splice_cigars)};" + + f"HAS_MODBASE_TAGS={int(has_modbase_tags)};" ) - logger.info(f"Checked (u)BAM headers in '{args.input_path}'.") + logger.info(f"Checked (u)BAM heads in '{args.input_path}'.") def argparser(): diff --git a/bin/workflow_glue/wfg_helpers/validators/wf.py b/bin/workflow_glue/wfg_helpers/validators/wf.py index 7975fab..d8903e5 100644 --- a/bin/workflow_glue/wfg_helpers/validators/wf.py +++ b/bin/workflow_glue/wfg_helpers/validators/wf.py @@ -5,7 +5,7 @@ import re from .default import SampleSheetValidator # noqa: ABS101 -_R_FORMULA_NAME_RE = re.compile(r"^[A-Za-z0-9][A-Za-z0-9._-]*$") +_R_FORMULA_NAME_RE = re.compile(r"^[A-Za-z][A-Za-z0-9_.]*$") def _split_covariates(value): @@ -31,8 +31,8 @@ def _validate_r_formula_names(names, label="Column"): details.append(f"Invalid names: {', '.join(invalid)}") raise ValueError( f"{label} names must be safe for R formulas. {'; '.join(details)}. " - "Names must start with a letter or number and contain " - "only letters, numbers, underscores, dots, and hyphens." + "Names must start with a letter and contain " + "only letters, numbers, underscores, and dots." ) diff --git a/bin/workflow_glue_r/R/common.R b/bin/workflow_glue_r/R/common.R index b41590d..4a87e8d 100644 --- a/bin/workflow_glue_r/R/common.R +++ b/bin/workflow_glue_r/R/common.R @@ -1,24 +1,3 @@ -workflow_glue_r_is_r_formula_name <- function(name) { - is.character(name) && - length(name) == 1 && - grepl("^[A-Za-z][A-Za-z0-9_.]*$", name) -} - -workflow_glue_r_validate_r_formula_names <- function(names, label = "Column") { - invalid <- names[!vapply(names, workflow_glue_r_is_r_formula_name, logical(1))] - if (length(invalid) > 0) { - stop( - sprintf( - "%s names must be safe for R formulas. Invalid names: %s. Names must start with a letter and contain only letters, numbers, underscores, and dots.", - label, - paste(invalid, collapse = ", ") - ), - call. = FALSE - ) - } - invisible(names) -} - workflow_glue_r_normalise_tsv_value <- function(value) { if (length(value) == 0 || all(is.na(value))) { return(NA_character_) diff --git a/bin/workflow_glue_r/R/de_analysis.R b/bin/workflow_glue_r/R/de_analysis.R index 80e5105..d62994d 100644 --- a/bin/workflow_glue_r/R/de_analysis.R +++ b/bin/workflow_glue_r/R/de_analysis.R @@ -59,10 +59,6 @@ de_analysis_arg_parser <- function() { de_validate_inputs <- function(tx_se, gene_se, sample_df, argv) { covariates <- workflow_glue_r_parse_csv_list(argv$covariates) - workflow_glue_r_validate_r_formula_names( - c(argv$condition_column, covariates), - label = "Design column" - ) if (any(duplicated(colnames(tx_se)))) { stop("Transcript RDS sample names must be unique.", call. = FALSE) diff --git a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R b/bin/workflow_glue_r/tests/testthat/test_de_analysis.R index a4c042d..80c9b09 100644 --- a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R +++ b/bin/workflow_glue_r/tests/testthat/test_de_analysis.R @@ -13,39 +13,6 @@ testthat::test_that("covariates parsed and trimmed", { ) }) -# R formulas break on spaces/hyphens in column names (e.g., ~ `time point` produces errors). -# Validate condition_column and covariate names are safe (alphanumeric + underscore). -testthat::test_that("formula-unsafe column names rejected", { - tx_se <- make_test_tx_se() - gene_se <- make_test_gene_se() - sample_df <- data.frame( - alias = colnames(tx_se), - condition = rep(c("control", "treated"), each = 3), - batch = rep(c("b1", "b2", "b1"), 2), - stringsAsFactors = FALSE - ) - argv <- list( - transcript_rds = "transcripts.rds", - gene_rds = "genes.rds", - sample_sheet = "sample_sheet.csv", - condition_column = "time point", - covariates = "batch", - reference_level = NULL - ) - argv <- workflow_glue_r_normalise_args(argv, de_analysis_arg_spec()) - - testthat::expect_error( - de_validate_inputs(tx_se, gene_se, sample_df, argv), - "Design column names must be safe for R formulas" - ) - - argv$condition_column <- "condition" - argv$covariates <- "batch-id" - testthat::expect_error( - de_validate_inputs(tx_se, gene_se, sample_df, argv), - "Invalid names: batch-id" - ) -}) testthat::test_that("DTU transcript output renames contrast-specific log2fold column", { contrast_name <- "treated_control" diff --git a/docs/04_install_and_run.md b/docs/04_install_and_run.md index f6995b4..4105bba 100644 --- a/docs/04_install_and_run.md +++ b/docs/04_install_and_run.md @@ -37,3 +37,21 @@ the following command: ``` nextflow pull epi2me-labs/wf-transcriptomes ``` + +A demo dataset is provided for testing of the workflow. +It can be downloaded and unpacked using the following commands: +``` +wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/v2/wf-transcriptomes-demo.tar.gz +tar -xzvf wf-transcriptomes-demo.tar.gz +``` +The workflow can then be run with the downloaded demo data using: +``` +nextflow run epi2me-labs/wf-transcriptomes \ + --bam 'wf-transcriptomes-demo/samples' \ + --de_analysis \ + --direct_rna \ + --ref_annotation 'wf-transcriptomes-demo/gencode.v22.annotation.chr20.gtf' \ + --ref_genome 'wf-transcriptomes-demo/hg38_chr20.fa' \ + --sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv' \ + -profile standard +``` diff --git a/docs/07_input_parameters.md b/docs/07_input_parameters.md index 75ac67d..0a86add 100644 --- a/docs/07_input_parameters.md +++ b/docs/07_input_parameters.md @@ -48,7 +48,6 @@ | 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`. | | | force_alignment | boolean | Force re-alignment of input BAM files. | Read alignment is skipped if the existing sequence names in the aligned BAM match the provided reference. Enable this if the existing alignments used incorrect minimap2 presets (e.g. missing --splice or direct RNA settings). | False | | ndr | number | Optional bambu novel discovery rate override. | Lower values are more conservative (higher precision), while higher values are more permissive (higher novel-discovery sensitivity). See the [`bambu` repository](https://github.com/GoekeLab/bambu) for method details. | | diff --git a/evaluation/tests.sh b/evaluation/tests.sh deleted file mode 100755 index f6b95a7..0000000 --- a/evaluation/tests.sh +++ /dev/null @@ -1,60 +0,0 @@ -#!/usr/bin/env bash - -# A few simple tests with different combinations of CLI options -# Run from within an appropriate active conda environment - -if [[ "$#" -lt 1 ]]; then - echo "usage: tests.sh [nextflow.config]" - exit 1 -fi - -if [[ "$#" -eq 1 ]]; then - config='' -fi - -if [[ "$#" -eq 2 ]]; then - config="-c $2"; -fi - -SCRIPT_DIR="$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )" -cd $SCRIPT_DIR/../; - -singledir="test_data/fastq" -multisampledir="test_data/demultiplexed_fastq" - -# This is for when using SIRV dataset with non-canonical spice junctions -#"--minimap2_opts '-uf --splice-flank=no'" -results=() - -# Reference based tests -OUTPUT=$1/reference_single_dir; -nextflow run . --fastq $singledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no' \ ---ref_annotation test_data/SIRV_isofroms.gtf -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace -resume; -r=$? -results+=("$(basename $OUTPUT): $r") - -OUTPUT=$1/multiple_samples; -nextflow run . --fastq $multisampledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no'\ ---ref_annotation test_data/SIRV_isofroms.gtf -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \ ---sample_sheet test_data/sample_sheet -resume; -r=$? -results+=("$(basename $OUTPUT): $r") - -OUTPUT=$1/reference_no_ref_annotation; -nextflow run . --fastq $singledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no'\ - -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace -resume; -r=$? -results+=("$(basename $OUTPUT): $r") - -# Force split_bam to make multiple alignment bundles -OUTPUT=$1/reference_frce_split_bam; -nextflow run . --fastq $singledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no'\ ---ref_annotation test_data/SIRV_isofroms.gtf -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \ ---bundle_min_reads 5 -resume; -r=$? -results+=("$(basename $OUTPUT): $r") - -echo "Exit status codes for each test" -for value in "${results[@]}"; do - echo "${value}" -done \ No newline at end of file diff --git a/lib/common.nf b/lib/common.nf index 1aed50d..aa2a9e8 100644 --- a/lib/common.nf +++ b/lib/common.nf @@ -67,7 +67,7 @@ process minimap2_alignment { 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'} + errorStrategy "retry" // retry on any error as we may not exit with appropriate codes when exceeding memory 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] diff --git a/lib/ingress.nf b/lib/ingress.nf index d457634..f10439f 100644 --- a/lib/ingress.nf +++ b/lib/ingress.nf @@ -340,6 +340,9 @@ def fastq_ingress(Map arguments, aln_ref_ch = null) `cram` outputs a CRAM file with index (.cram, .crai) * - "sort_threads": number of threads to use for sorting aligned reads (default: 1) * - "force_alignment": boolean. Run alignment regardless of inputs. + * - "checked_bam_hook": closure. Optional hook called after all BAM header/read checks + * complete. Receives a list of `[meta, paths]` tuples and must return replacement + * tuples of the same shape. Can warn, error, or add routing metadata. * @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]`. @@ -361,8 +364,13 @@ def xam_ingress(Map arguments, aln_ref_ch = null) "fastcat_extra_args": "", "fastq_chunk": null, "force_alignment": false, + "checked_bam_hook": null, ] ) + // immediately separate hooks from margs as they cannot be serialised + // and are currently only needed in this scope + def checked_bam_hook = margs.remove("checked_bam_hook") + margs["fastq_chunk"] ?= 0 // cant pass null through channel if (margs["return_fastq"] && aln_ref_ch){ @@ -408,14 +416,24 @@ def xam_ingress(Map arguments, aln_ref_ch = null) [meta + [src_xam: src_xam, src_xai: src_xai], paths] } | combine(aln_ref) - | checkBamHeaders - | map { meta, paths, is_unaligned_env, mixed_sq_headers_env, is_sorted_env, has_reads_env -> + | checkBamHeads + | map { + meta, + paths, + is_unaligned_env, + mixed_sq_headers_env, + is_sorted_env, + has_reads_env, + has_splice_cigars_env, + has_modbase_tags_env -> // convert the env. variables from strings ('0' or '1') into bools boolean is_unaligned = is_unaligned_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 + boolean has_splice_cigars = has_splice_cigars_env as int as boolean + boolean has_modbase_tags = has_modbase_tags_env as int as boolean // throw an error if there was a sample with mixed headers if (mixed_sq_headers) { error "Found mixed headers in (u)BAM files of sample '${meta.alias}'." @@ -423,7 +441,24 @@ def xam_ingress(Map arguments, aln_ref_ch = null) // 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, has_reads: has_reads], paths] + [ + meta + [ + is_unaligned: is_unaligned, + is_sorted: is_sorted, + has_reads: has_reads, + has_splice_cigars: has_splice_cigars, + has_modbase_tags: has_modbase_tags, + ], + paths + ] + } + + if (aln_ref_ch && checked_bam_hook) { + ch_check_bams = ch_check_bams + | collect(flat: false) + | flatMap { checked_bams -> + checked_bam_hook(checked_bams) + } } // Handle alignment @@ -431,7 +466,7 @@ def xam_ingress(Map arguments, aln_ref_ch = null) alignment_fork = ch_check_bams | branch { meta, paths -> - to_align: (meta.is_unaligned || margs.force_alignment) && meta.has_reads + to_align: (meta.is_unaligned || meta.requires_alignment || margs.force_alignment) && meta.has_reads noalign: true } if (margs["output_xam_fmt"] == "bam"){ @@ -448,14 +483,18 @@ def xam_ingress(Map arguments, aln_ref_ch = null) // 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' } + // remove alignment routing metadata that is no longer required + def newmeta = meta.findAll { + k, v -> !(k in ['has_reads', 'requires_alignment']) + } [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' } + def newmeta = meta.findAll { + k, v -> !(k in ['has_reads', 'requires_alignment']) + } [newmeta, paths] } } else { @@ -769,7 +808,7 @@ process fastcat { """ } -process checkBamHeaders { +process checkBamHeads { label "ingress" label "wf_common" cpus 1 @@ -784,6 +823,8 @@ process checkBamHeaders { env(MIXED_SQ_HEADERS), env(IS_SORTED), env(HAS_READS), + env(HAS_SPLICE_CIGARS), + env(HAS_MODBASE_TAGS), ) script: String ref_arg = ref.fileName.name == OPTIONAL_FILE.name ? "" : "--ref $ref --ref_idx $ref_index" diff --git a/main.nf b/main.nf index 7dabfd7..4824c89 100644 --- a/main.nf +++ b/main.nf @@ -20,6 +20,8 @@ process getVersions { publishDir "${params.out_dir}", mode: 'copy', pattern: "versions.txt" cpus 1 memory "2 GB" + input: + path "additional_versions.txt" output: path "versions.txt" script: @@ -28,6 +30,36 @@ process getVersions { samtools --version | head -n 1 | sed 's/ /,/' >> versions.txt gffread --version | sed 's/^/gffread,/' >> versions.txt || true Rscript -e 'pkgs <- c("bambu", "DESeq2", "DEXSeq"); for (pkg in pkgs) {cat(pkg, as.character(packageVersion(pkg)), sep = ","); cat("\\n")}' >> versions.txt + cat additional_versions.txt >> versions.txt + """ +} + + +process getSqantiVersion { + label "wf_transcriptomes_sqanti" + cpus 1 + memory "2 GB" + output: + path "versions.txt" + script: + """ + sqanti3_qc.py --version | sed 's/ /,/' >> versions.txt + """ +} + + +process getModkitVersion { + label "modkit" + cpus 1 + memory "2 GB" + input: + path "old_versions.txt" + output: + path "versions.txt" + script: + """ + cp old_versions.txt versions.txt + modkit --version | sed 's/ /,/' >> versions.txt """ } @@ -100,7 +132,7 @@ workflow wf { ref_genome ref_annotation main: - software_versions = getVersions() + software_versions = getVersions(getModkitVersion(getSqantiVersion())) workflow_params = getParams() transcriptome_results = transcriptome(reads, ref_genome, ref_annotation, sample_sheet) @@ -235,6 +267,70 @@ workflow { "sort_threads": 3, ] + def checked_bam_hook = { checked_bams -> + if (params.force_alignment) { + // short-circuit this hook to allow for forced realignment + return checked_bams + } + + // determine BAMs with and without splice-aware CIGAR evidence + ArrayList aliases_with_splice_cigars = [] + ArrayList aliases_without_splice_cigars = [] + checked_bams.each { + def meta = it[0] + if (meta.has_reads) { + if (meta.has_splice_cigars) { + aliases_with_splice_cigars.add(meta.alias) + } else { + aliases_without_splice_cigars.add(meta.alias) + } + } + } + boolean any_with_splice_cigars = aliases_with_splice_cigars + boolean any_without_splice_cigars = aliases_without_splice_cigars + + // explode on mixed splice-aware status + if (any_with_splice_cigars && any_without_splice_cigars) { + log.error """ + Cannot proceed with mixed splice-aware CIGAR evidence in input BAMs. + + This workflow requires splice-aware alignment. It cannot safely realign + only the BAMs that appear to lack splice-aware CIGARs, because the + original alignment parameters for the other BAMs are unknown. + + It is unlikely for a splice-aware aligned transcriptome BAM to have no + reads with N CIGAR operations. Inspect your BAMs to understand why their + splice-aware alignment evidence differs. Either remove the BAMs without + splice-aware CIGARs from the analysis, or use --force_alignment to realign + all inputs. + + Samples with splice CIGARs: ${aliases_with_splice_cigars.join(', ')} + Samples without splice CIGARs: ${aliases_without_splice_cigars.join(', ')} + """.stripIndent().trim() + error "Mixed splice-aware CIGAR evidence in input BAMs." + } + + // all read-containing BAMs have splice-aware CIGAR evidence, or there are no reads + if (!any_without_splice_cigars) { + return checked_bams + } + + // if we're here, all read-containing BAMs lack splice-aware CIGAR evidence + // fiddle the metamap to set requires_alignment (if there are reads to align) + log.warn "No input BAMs appear to contain splice-aware CIGAR evidence. " + + "The workflow will realign all inputs." + checked_bams.collect { + def meta = it[0] + def paths = it[1] + [ + meta + [ + requires_alignment: meta.has_reads, + ], + paths + ] + } + } + if (params.fastq) { samples = fastq_ingress([ "input": params.fastq, @@ -243,6 +339,7 @@ workflow { samples = xam_ingress([ "input": params.bam, "force_alignment": params.force_alignment, + "checked_bam_hook": checked_bam_hook, ] + ingress_args, ref_genome) } diff --git a/nextflow.config b/nextflow.config index a1cdbb8..5bfea6a 100644 --- a/nextflow.config +++ b/nextflow.config @@ -22,7 +22,6 @@ params { analyse_unclassified = false analyse_fail = false - threads = 4 mod_codes = null force_alignment = false ndr = null @@ -44,9 +43,9 @@ params { per_read_stats = false allow_multiple_basecall_models = false example_cmd = [ + "--bam 'wf-transcriptomes-demo/samples'", "--de_analysis", "--direct_rna", - "--fastq 'wf-transcriptomes-demo/differential_expression_fastq'", "--ref_annotation 'wf-transcriptomes-demo/gencode.v22.annotation.chr20.gtf'", "--ref_genome 'wf-transcriptomes-demo/hg38_chr20.fa'", "--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv'", @@ -68,7 +67,7 @@ manifest { description = 'Long-read transcript discovery, quantification, differential expression, QC and mod counting.' mainScript = 'main.nf' nextflowVersion = '>=23.04.2' - version = 'v2.0.0-rc1' + version = 'v2.0.0' } process { diff --git a/nextflow_schema.json b/nextflow_schema.json index 1bb5adf..e864d2f 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -4,6 +4,8 @@ "title": "epi2me-labs/wf-transcriptomes", "workflow_title": "Transcriptomes", "description": "Long-read transcript discovery, quantification, differential expression, QC and mod counting.", + "demo_url": "https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/v2/wf-transcriptomes-demo.tar.gz", + "aws_demo_url": "https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/wf-transcriptomes-demo/transcriptomes.v2.aws.nextflow.config", "url": "https://github.com/epi2me-labs/wf-transcriptomes", "type": "object", "definitions": { @@ -167,11 +169,6 @@ "type": "object", "description": "Performance tuning and tool-specific advanced parameters.", "properties": { - "threads": { - "type": "integer", - "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.", diff --git a/subworkflows/differential_expression.nf b/subworkflows/differential_expression.nf index 645fc4c..f7bbb90 100644 --- a/subworkflows/differential_expression.nf +++ b/subworkflows/differential_expression.nf @@ -3,7 +3,7 @@ nextflow.enable.dsl = 2 process runDifferentialAnalysis { label "wf_transcriptomes" - cpus { params.threads ?: 4 } + cpus 4 memory "32 GB" input: path transcript_rds diff --git a/subworkflows/mods.nf b/subworkflows/mods.nf index e00247e..ebd8f77 100644 --- a/subworkflows/mods.nf +++ b/subworkflows/mods.nf @@ -1,35 +1,6 @@ -// 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 } + cpus 4 memory "16 GB" input: tuple val(alias), @@ -158,15 +129,12 @@ workflow mods { 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} + xams.branch { + nomods: it[0].has_modbase_tags != true + return it[0].alias + mods: it[0].has_modbase_tags == true + return [it[0].alias] + it // prepend alias for joining + }.set{xams_with} // warn for samples without mods xams_with.nomods.subscribe { diff --git a/subworkflows/transcriptome.nf b/subworkflows/transcriptome.nf index 846c476..41f5c20 100644 --- a/subworkflows/transcriptome.nf +++ b/subworkflows/transcriptome.nf @@ -138,7 +138,7 @@ process buildSampleTranscriptomeFasta { process runJointSqanti { label "wf_transcriptomes_sqanti" - cpus { params.threads ?: 4 } + cpus 2 memory "24 GB" publishDir "${params.out_dir}/cohort", mode: "copy", @@ -172,7 +172,7 @@ process runJointSqanti { process runPerSampleSqanti { label "wf_transcriptomes_sqanti" - cpus { params.threads ?: 4 } + cpus 2 memory "24 GB" publishDir "${params.out_dir}/samples/${meta.alias}", mode: "copy",