Merge branch 'v200-rc2' into 'dev'
wf-transcriptomes v2.0.0 See merge request epi2melabs/workflows/wf-transcriptomes!320
This commit is contained in:
Sam Nicholls
commit
1fa032ae0e
@ -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:
|
||||
|
||||
@ -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.
|
||||
|
||||
19
README.md
19
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. | |
|
||||
|
||||
@ -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
|
||||
@ -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.")
|
||||
|
||||
|
||||
@ -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)",
|
||||
]
|
||||
|
||||
|
||||
|
||||
@ -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"))
|
||||
@ -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
|
||||
|
||||
@ -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():
|
||||
|
||||
@ -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."
|
||||
)
|
||||
|
||||
|
||||
|
||||
@ -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_)
|
||||
|
||||
@ -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)
|
||||
|
||||
@ -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"
|
||||
|
||||
@ -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
|
||||
```
|
||||
|
||||
@ -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. | |
|
||||
|
||||
@ -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 <outdir> [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
|
||||
@ -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]
|
||||
|
||||
@ -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"
|
||||
|
||||
99
main.nf
99
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)
|
||||
}
|
||||
|
||||
|
||||
@ -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 {
|
||||
|
||||
@ -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.",
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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 {
|
||||
|
||||
@ -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",
|
||||
|
||||
Loading…
Reference in New Issue
Block a user