kfukuta/wf-transcriptomes-v202
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'template-ingress' into 'dev'

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

See merge request epi2melabs/workflows/wf-transcriptomes!247
This commit is contained in:
Sam Nicholls 2026-05-14 09:12:38 +00:00
commit 42d29d3a4d
23 changed files with 942 additions and 658 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
.gitlab-ci.yml
View File

@ -15,6 +15,7 @@ variables:
PYTEST_CONTAINER_CONFIG_KEY: "common_sha"
RTEST_CONTAINER_NAME: "wf-transcriptomes-core"
RTEST_CONTAINER_CONFIG_KEY: "container_sha"
WF_TEMPLATE_ENFORCEMENT_BRANCH: "CW-6552"
macos-run:
# Let's avoid those ARM64 runners for now
@ -56,7 +57,7 @@ docker-run:
parallel:
matrix:
- MATRIX_NAME: [
"de", "discover", "igv",
"de-poscounts-fallback", "discover", "igv",
"smoke_discover", "smoke_fixed", "smoke_direct_rna", "smoke_de",
"no_annotation", "invalid_mode", "conflicting_flags"
]
@ -70,18 +71,17 @@ docker-run:
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq \
--ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf --cdna_preprocess --pychopper_backend phmm"
NF_IGNORE_PROCESSES: preprocess_reads,filter_unstranded_annotation,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
--ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf"
NF_IGNORE_PROCESSES: filter_unstranded_annotation,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "no_annotation"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq \
--ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: preprocess_reads
ASSERT_NEXTFLOW_FAILURE: "1"
ASSERT_NEXTFLOW_FAILURE_REXP: "Missing required parameter: --ref_annotation"
- if: $MATRIX_NAME == "de"
- if: $MATRIX_NAME == "de-poscounts-fallback"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
@ -90,7 +90,9 @@ docker-run:
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gtf \
--direct_rna --minimap2_index_opts '-k 15' --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression/sample_sheet.csv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: preprocess_reads,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
NF_IGNORE_PROCESSES: faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
AFTER_NEXTFLOW_CMD: >
grep -Eq '"deseq2_size_factor_method": "poscounts"' ${CI_PROJECT_NAME}/de_analysis/de_qc_stats.json
- if: $MATRIX_NAME == "only_differential_expression"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
@ -103,13 +105,12 @@ docker-run:
--igv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: >
preprocess_reads,gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,decompress_ref,
gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,decompress_ref,
build_minimap_index,validate_ref_annotation,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "smoke_discover"
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf"
NF_IGNORE_PROCESSES: preprocess_reads
AFTER_NEXTFLOW_CMD: >
test -f ${CI_PROJECT_NAME}/cohort/transcripts.gtf &&
test -f ${CI_PROJECT_NAME}/cohort/cohort.transcriptome.fa &&
@ -118,7 +119,6 @@ docker-run:
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --transcriptome_mode fixed_annotation"
NF_IGNORE_PROCESSES: preprocess_reads
AFTER_NEXTFLOW_CMD: >
test -f ${CI_PROJECT_NAME}/cohort/transcripts.gtf &&
test -f ${CI_PROJECT_NAME}/cohort/transcript_counts.tsv
@ -126,15 +126,13 @@ docker-run:
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --direct_rna"
NF_IGNORE_PROCESSES: preprocess_reads
AFTER_NEXTFLOW_CMD: >
test -f ${CI_PROJECT_NAME}/cohort/alignments/sampleA.aligned.sorted.bam &&
test -f ${CI_PROJECT_NAME}/cohort/alignments/sampleA/reads.bam &&
test -f ${CI_PROJECT_NAME}/cohort/sqanti_cohort/classification_summary.tsv
- if: $MATRIX_NAME == "smoke_de"
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/de --sample_sheet test_data/smoke/sample_sheet_de.csv --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --de_analysis --reference_level control --covariates batch"
NF_IGNORE_PROCESSES: preprocess_reads
AFTER_NEXTFLOW_CMD: >
test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dge.tsv &&
test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dtu_transcript.tsv
@ -142,12 +140,5 @@ docker-run:
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --transcriptome_mode nonsense"
NF_IGNORE_PROCESSES: preprocess_reads
ASSERT_NEXTFLOW_FAILURE: "1"
ASSERT_NEXTFLOW_FAILURE_REXP: "nonsense is not a valid choice"
- if: $MATRIX_NAME == "conflicting_flags"
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --direct_rna --cdna_preprocess"
ASSERT_NEXTFLOW_FAILURE: "1"
ASSERT_NEXTFLOW_FAILURE_REXP: ".*--cdna_preprocess cannot be used together with --direct_rna"

57
README.md
View File

@ -18,7 +18,6 @@ analysis.
The workflow supports:
+ transcript identification from either cDNA or direct RNA reads
+ optional cDNA preprocessing with `pychopper` before alignment
+ transcript discovery guided by a supplied genome and annotation
+ quantification against a supplied reference annotation
+ optional transcript classification and QC with `SQANTI3`
@ -211,15 +210,7 @@ barcode04,treated_rep2,test_sample,treated,b2
This example is suitable for a multiplexed run and also satisfies the minimum
requirements for a two-group DE/DTU comparison.
### 3. Optional cDNA preprocessing
When `--cdna_preprocess` is enabled for cDNA libraries, the workflow runs
`pychopper` before alignment to classify, orient, and trim full-length reads.
This preprocessing stage is controlled by `--cdna_kit`,
`--pychopper_backend`, and optional extra `--pychopper_opts`, and its outputs
are published alongside the ingress results for each sample.
### 4. Genome alignment
### 3. Genome alignment
Each sample is aligned to the supplied reference genome with
[`minimap2`](https://github.com/lh3/minimap2), then sorted and indexed with
@ -227,14 +218,14 @@ Each sample is aligned to the supplied reference genome with
`cohort/alignments/` are the main alignment files used for transcriptome
analysis, optional `SQANTI3` QC, and optional IGV viewing.
### 5. Cohort transcriptome construction
### 4. Cohort transcriptome construction
All aligned samples are analysed together with `bambu` to produce the primary
cohort transcriptome, transcript counts, gene counts, and the `RDS` objects used
for downstream differential analysis. This shared model is the main cohort-level
result and is published under `cohort/`.
### 6. Independent per-sample transcriptomes
### 5. Independent per-sample transcriptomes
Each sample is also processed separately with `bambu` so the workflow produces
sample-specific GTF, FASTA, count tables, and metadata under
@ -242,7 +233,7 @@ sample-specific GTF, FASTA, count tables, and metadata under
specific transcript models without changing the shared cohort transcriptome used
for DE/DTU.
### 7. Transcript sequence generation and QC
### 6. Transcript sequence generation and QC
Transcript FASTA files are derived from GTF plus genome using `gffread`.
When `--skip_sqanti` is not set, `SQANTI3` classifies the cohort and per-sample
@ -250,7 +241,7 @@ transcriptomes and produces structural QC summaries. The cohort `SQANTI3`
results live under `cohort/sqanti_cohort/`, while per-sample `SQANTI3`
directories are published under `samples/<alias>/<alias>_sqanti/`.
### 8. Optional DE and DTU analysis
### 7. Optional DE and DTU analysis
When `--de_analysis` is enabled, the workflow checks the experimental design,
runs `DESeq2` for differential gene expression, and runs `DEXSeq` for
@ -258,7 +249,7 @@ differential transcript usage. These analyses use the shared `bambu` outputs
and the design columns in the sample sheet, and each comparison is written to
its own subdirectory under `de_analysis/<contrast>/`.
### 9. What you need to provide
### 8. What you need to provide
The workflow's analysis is controlled by a user provided genome, annotation, and
`bambu` mode.
@ -270,13 +261,10 @@ The workflow's analysis is controlled by a user provided genome, annotation, and
* `--ref_transcriptome` has been removed; if you want annotation-based
quantification, use `--transcriptome_mode fixed_annotation` together with
`--ref_genome` and `--ref_annotation`
* `--cdna_preprocess` enables the `pychopper` cDNA preprocessing stage, and
`--cdna_kit`, `--pychopper_backend`, and `--pychopper_opts` control that stage
* when `--de_analysis` is enabled, the sample sheet must contain `alias`, the
primary condition column, and any requested columns named in `--covariates`
* `--cdna_preprocess` must not be combined with `--direct_rna`
### 10. How to read the output folder
### 9. How to read the output folder
The published outputs are organised around a small number of top-level
directories:
@ -315,7 +303,6 @@ directories:
| ref_annotation | string | Reference transcript annotation in GTF or GFF format. | Required in both discover and fixed_annotation modes. | |
| transcriptome_mode | string | How bambu should prepare the transcriptome model. | Use discover for reference-guided transcript discovery and quantification, or fixed_annotation for quantification only against the supplied annotation. | discover |
| direct_rna | boolean | Set this for direct RNA sequencing libraries. | | False |
| cdna_preprocess | boolean | Apply the optional cDNA preprocessing stage before alignment. | This is only relevant for cDNA libraries and must not be combined with `--direct_rna`. | False |
### Sample Options
@ -351,9 +338,6 @@ directories:
| threads | integer | Thread count to use for the core workflow processes. | | 4 |
| minimap2_opts | string | Extra command-line options to pass to minimap2. | | |
| ndr | number | Optional bambu novel discovery rate override. | | |
| cdna_kit | string | ONT cDNA kit identifier used for pychopper preprocessing. | Only used when `--cdna_preprocess` is enabled. The workflow derives the pychopper kit code from this value. | SQK-PCS109 |
| pychopper_backend | string | Primer-detection backend to use for pychopper preprocessing. | | edlib |
| pychopper_opts | string | Extra command-line options to pass to pychopper. | | |
| skip_sqanti | boolean | Skip SQANTI3 transcript classification and QC. | | False |
| sqanti_skip_orf | boolean | Skip ORF prediction during SQANTI3 QC. | | True |
| sqanti_extra_args | string | Extra command-line options to pass to SQANTI3. | | |
@ -374,9 +358,9 @@ Output files may be aggregated including information for all samples or provided
| Per-read stats | ingress_results/{{ alias }}/fastcat_stats/per-read-stats.tsv.gz | Read statistics for individual reads in a sample, when this output is enabled. | per-sample |
| Ingress reads | ingress_results/{{ alias }}/seqs.fastq.gz | Reads prepared from the input data for downstream analysis. | per-sample |
| Ingress metadata | ingress_results/{{ alias }}/metamap.json | Per-sample metadata used by the workflow. | per-sample |
| Aligned BAM | cohort/alignments/{{ alias }}.aligned.sorted.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
| Aligned BAM index | cohort/alignments/{{ alias }}.aligned.sorted.bam.bai | Index for the aligned BAM. | per-sample |
| Alignment summary | cohort/alignments/{{ alias }}.flagstat.txt | samtools flagstat output for the aligned BAM. | per-sample |
| Aligned BAM | cohort/alignments/{{ alias }}/reads.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
| Aligned BAM index | cohort/alignments/{{ alias }}/reads.bam.bai | Index for the aligned BAM. | per-sample |
| Alignment summary | cohort/alignments/{{ alias }}/bamstats.flagstat.tsv | bamstats flagstat summary for the aligned BAM. | per-sample |
| Reference and annotation preparation summary | cohort/reference/annotation_reference_summary.json | Summary of reference and annotation preparation, including seqname overlap, build/provider hints, and excluded unstranded annotation counts. | aggregated |
| Excluded unstranded annotation records | cohort/reference/unstranded_annotation.gtf | Full set of annotation records excluded because their strand was not '+' or '-'. Present only when unstranded records are found. | aggregated |
| Cohort transcriptome GTF | cohort/transcripts.gtf | Joint bambu transcript model used as the primary cohort transcriptome. | aggregated |
@ -429,7 +413,6 @@ Find related RNA and cDNA sequencing protocols in the
primary condition column, and any columns named in `--covariates`.
+ DE/DTU requires at least two condition levels and at least two samples per
level.
+ `--cdna_preprocess` must not be combined with `--direct_rna`.
+ See how to interpret common Nextflow exit codes
[here](https://labs.epi2me.io/trouble-shooting/).
@ -464,17 +447,6 @@ setting that chose how the workflow behaved. In the current version,
The options `--transcriptome_mode discover` or `--transcriptome_mode
fixed_annotation` should be used to choose between the modes of operation.
#### I expected pychopper-style preprocessing for cDNA libraries
cDNA preprocessing is available through
`--cdna_preprocess`, which runs `pychopper` before alignment. The
parameters for that stage are `--cdna_kit`, `--pychopper_backend`, and
`--pychopper_opts`.
Enable `--cdna_preprocess` for cDNA libraries when you want
`pychopper` preprocessing, or use `--direct_rna` without
`--cdna_preprocess` for direct RNA data.
#### I cannot find the old flat DE output files
In the previous workflow version, DE files appeared directly under `de_analysis/`.
@ -491,9 +463,7 @@ In the current versions, the output folder is organised around `ingress_results/
`igv_reference/`.
Primary shared transcriptome results are under `cohort/`, and
`samples/<alias>/` for sample-specific models. If `--cdna_preprocess` is enabled,
look under `ingress_results/<alias>/` for the `pychopper` outputs associated
with each sample.
`samples/<alias>/` for sample-specific models.
#### My DE/DTU run fails because of `--sample_sheet`
@ -507,11 +477,6 @@ are required, and each level must contain at least two samples.
What to change: verify the sample sheet columns first, then check
`--condition_column`, `--covariates`, and `--reference_level`.
#### I combined `--direct_rna` with `--cdna_preprocess`
This is not a valid combination of parameters. Use `--direct_rna` for direct RNA
libraries and leave `--cdna_preprocess` unset.
#### I hit genome/annotation validation or strand-related annotation warnings
The workflow validates that the annotation and genome share

6
bin/workflow_glue/models/common.py
View File

@ -206,10 +206,11 @@ class Sample:
self.additional_identifiers.append(sample_identifier)
return self.additional_identifiers
# TODO: Integrate this into WorkflowBaseModel
def to_json(self, filename):
"""Save class as JSON."""
with open(filename, 'w') as f:
json.dump(asdict(self), f, default=str, indent=2)
json.dump(asdict(self), f, default=str, indent=2, allow_nan=False)
def get_reportable_qc_status(self, max_criteria=4):
"""Store global status of the sample and list of QC criteria to show.
@ -364,7 +365,8 @@ class WorkflowResult(WorkflowBaseModel):
self.versions = versions
return self.versions
# TODO: Integrate this into WorkflowBaseModel
def to_json(self, filename):
"""Save class as JSON."""
with open(filename, 'w') as f:
json.dump(asdict(self), f, default=str, indent=2)
json.dump(asdict(self), f, default=str, indent=2, allow_nan=False)

41
bin/workflow_glue/report.py
View File

@ -81,17 +81,6 @@ def _sqanti_tables(sqanti_dir):
return tables
def _pychopper_tables(pychopper_dir):
tables = {}
for summary in sorted(Path(pychopper_dir).rglob("pychopper_summary.tsv")):
table = _read_table(summary)
if table is None or table.empty:
continue
label = summary.parent.name.replace("_pychopper_output", "")
tables[label] = table
return tables
def _top_results(de_dir, filename, n=20):
tables = {}
for contrast_dir in sorted(Path(de_dir).iterdir()):
@ -193,6 +182,7 @@ def _collect_de_method_rows(de_qc):
fallback = contrast_data.get("deseq2_dispersion_fallback") or {}
fallback_applied = bool(fallback.get("applied", False))
deseq2_method = fallback.get("method_used")
deseq2_size_factors = contrast_data.get("deseq2_size_factor_method") or "ratio"
if not deseq2_method:
deseq2_method = "gene-wise" if fallback_applied else "parametric"
@ -201,6 +191,7 @@ def _collect_de_method_rows(de_qc):
deseq2_gene_wise.append(contrast_name)
dexseq_method = contrast_data.get("dexseq_dispersion_method") or "parametric"
dexseq_size_factors = contrast_data.get("dexseq_size_factor_method") or "ratio"
if dexseq_method == "gene-wise":
dexseq_gene_wise.append(contrast_name)
@ -213,11 +204,13 @@ def _collect_de_method_rows(de_qc):
rows.append(
{
"Contrast": contrast_name,
"DESeq2 size factors": deseq2_size_factors,
"DESeq2 dispersion": (
f"{deseq2_method} (fallback)"
if fallback_applied
else deseq2_method
),
"DEXSeq size factors": dexseq_size_factors,
"DEXSeq dispersion": dexseq_method,
"DEXSeq covariates dropped": (
", ".join(dropped_covariates) if dropped_covariates else "none"
@ -257,10 +250,19 @@ def main(args):
sample_names = tuple(
item["alias"] for item in metadata if item.get("has_stats")
)
flagstats = tuple(
Path(stats_dir) / "bamstats.flagstat.tsv" for stats_dir in stats
)
if len(stats) == 1:
stats = stats[0]
flagstats = flagstats[0]
sample_names = sample_names[0] if sample_names else None
fastcat.SeqSummary(stats, sample_names=sample_names)
fastcat.SeqSummary(
stats,
flagstat=flagstats,
sample_names=sample_names,
alignment_stats=True,
)
with report.add_section("Sample metadata", "Samples"):
tabs = Tabs()
@ -499,16 +501,6 @@ def main(args):
with tabs.add_tab(stats_file.stem.replace(".flagstat", "")):
pre(stats_file.read_text())
pychopper_tables = {}
if args.pychopper_dir and Path(args.pychopper_dir).exists():
pychopper_tables = _pychopper_tables(args.pychopper_dir)
if pychopper_tables:
with report.add_section("Pychopper preprocessing", "Pychopper"):
tabs = Tabs()
for label, table in pychopper_tables.items():
with tabs.add_tab(label):
DataTable.from_pandas(table, use_index=False)
sqanti_tables = _sqanti_tables(args.sqanti_dir)
if sqanti_tables:
with report.add_section("SQANTI3 classification", "SQANTI3"):
@ -853,11 +845,6 @@ def argparser():
required=True,
help="Per-sample output directory.",
)
parser.add_argument(
"--pychopper_dir",
default=None,
help="Pychopper output directory.",
)
parser.add_argument(
"--sqanti_dir",
required=True,

19
bin/workflow_glue/tests/test_report.py
View File

@ -217,25 +217,6 @@ def test_report_main_accepts_optional_file_sentinels(monkeypatch, tmp_path):
assert any("GRCh38" in table.to_string() for table in tables)
def test_pychopper_tables_uses_sample_directory_names(tmp_path):
"""Pychopper summaries should be keyed by sample alias."""
pychopper_dir = tmp_path / "pychopper"
sample_dir = pychopper_dir / "sampleA_pychopper_output"
sample_dir.mkdir(parents=True)
_write(
sample_dir / "pychopper_summary.tsv",
"Classification\tValue\nFull length\t10\nUnclassified\t2\n",
)
tables = report._pychopper_tables(pychopper_dir)
assert list(tables) == ["sampleA"]
assert list(tables["sampleA"]["Classification"]) == [
"Full length",
"Unclassified",
]
def test_report_main_renders_statistical_methods_and_warnings(
monkeypatch,
tmp_path,

151
bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py
View File

@ -1,5 +1,5 @@
"""Check (u)BAM files for `@SQ` lines whether they are the same in all headers."""
import os
from pathlib import Path
import sys
@ -8,58 +8,156 @@ import pysam
from ..util import get_named_logger, wf_parser # noqa: ABS101
def extract_header_info(xam_file, has_reads, check_ref=False):
"""Extract header information from a BAM/CRAM file."""
try:
f = pysam.AlignmentFile(xam_file, check_sq=False)
except (ValueError, IOError):
return None, None, None, False
with f:
# Extract SQ lines, comparing only SN/LN/M5 elements
# (see CW-4842 - ignore different SQ.UR values)
sq_lines = [{
"SN": sq["SN"],
"LN": sq["LN"],
"M5": sq.get("M5"),
} for sq in f.header.get("SQ", [])]
hd_lines = f.header.get("HD")
xam_reflen = None
if check_ref:
xam_reflen = set(zip(f.references, f.lengths))
# Whilst we have file open check for atleast 1 reads
# unless reads already found
if not has_reads:
try:
next(f.fetch(until_eof=True))
has_reads = True
except StopIteration:
has_reads = False
return sq_lines, hd_lines, xam_reflen, has_reads
def compare_ref_lengths(xam_reflen, ref_reflen, logger):
"""Compare reference lengths from FASTA and BAM/CRAM SQ lines."""
diff = ref_reflen.symmetric_difference(xam_reflen)
if len(diff) > 0:
rows = ["sequence_name sequence_length in_ref in_xam"]
rows += [
f"{reflen[0]} {reflen[1]} "
f"{'1' if reflen in ref_reflen else '0'} "
f"{'1' if reflen in xam_reflen else '0'}"
for reflen in diff
]
logger.info("\n".join(rows) + "\n")
logger.info(
" [WARN] In the input XAM there is at least one (name, length)"
" pair that does not map 1:1 between the input alignment"
" and reference. XAM will be realigned.\n"
)
return False
else:
logger.info(
"[OKAY] Input alignment and reference sequences map 1:1\n"
)
return True
def main(args):
"""Run the entry point."""
logger = get_named_logger("checkBamHdr")
# Check the reference if given.
ref_reflen = None
if args.ref:
try:
ref = pysam.FastaFile(
args.ref,
filepath_index=args.ref_idx)
except (ValueError, IOError):
logger.info(
f"[FAIL] {args.ref} reference input file could not be"
" read. Is it in the right format?"
)
sys.exit(os.EX_NOINPUT)
# get ref length for comparison with BAM/CRAM SQ lines
ref_reflen = set(zip(ref.references, ref.lengths))
if not args.input_path.is_dir():
raise ValueError(f"Input path '{args.input_path}' must be a directory.")
target_files = list(args.input_path.glob("*"))
target_files = sorted(list(args.input_path.glob("*")))
if not target_files:
raise ValueError(f"No files found in input directory '{args.input_path}'.")
# Loop over target files and check if there are `@SQ` lines in all headers or not.
# Set `is_unaligned` accordingly. If there are mixed headers (either with some files
# containing `@SQ` lines and some not or with different files containing different
# `@SQ` lines), set `mixed_headers` to `True`.
# Set `is_unaligned` accordingly.
# Detect mixed headers by comparing @SQ lines between files.
# Set mixed_sq_headers to True if they differ across files
# (either some files lack @SQ lines, or files have different @SQ lines).
# Also check if there is the SO line, to validate whether the file is (un)sorted.
first_sq_lines = None
mixed_headers = False
mixed_sq_headers = False
sorted_xam = False
xam_reflen = None
requires_realign = False
any_sq_lines = False # Track if any file has @SQ lines
has_reads = False
for xam_file in target_files:
# get the `@SQ` and `@HD` lines in the header
with pysam.AlignmentFile(xam_file, check_sq=False) as f:
# compare only the SN/LN/M5 elements of SQ to avoid labelling XAM with
# same reference but different SQ.UR as mixed_header (see CW-4842)
sq_lines = [{
"SN": sq["SN"],
"LN": sq["LN"],
"M5": sq.get("M5"),
} for sq in f.header.get("SQ", [])]
hd_lines = f.header.get("HD")
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:
# File couldn't be opened
logger.error(f"Failed to open {xam_file}")
continue
if sq_lines:
any_sq_lines = True
# Check if it is sorted.
# When there is more than one BAM, merging/sorting
# will happen regardless of this flag.
if hd_lines is not None and hd_lines.get('SO') == 'coordinate':
sorted_xam = True
if first_sq_lines is None:
# this is the first file
first_sq_lines = sq_lines
else:
# this is a subsequent file; check with the first `@SQ` lines
if sq_lines != first_sq_lines:
mixed_headers = True
mixed_sq_headers = True
break
# we set `is_unaligned` to `True` if there were no mixed headers and the last file
# didn't have `@SQ` lines (as we can then be sure that none of the files did)
is_unaligned = not mixed_headers and not sq_lines
# write `is_unaligned` and `mixed_headers` out so that they can be set as env.
# variables
if ref_reflen and not mixed_sq_headers and not requires_realign:
if sq_lines:
# if a reference was given and there are sq lines (indicating aligned)
# compare ref and xam header lengths to check alignment is to this ref
requires_realign = not compare_ref_lengths(
xam_reflen, ref_reflen, logger)
else:
requires_realign = True
# Assumption: presence of @SQ lines indicates aligned data.
# Set `is_unaligned` to `True` if there were no mixed headers and none of the
# files had `@SQ` lines, or if found to be incorrectly aligned based on
# SQ/ref length comparison.
is_unaligned = (not mixed_sq_headers and not any_sq_lines) or requires_realign
# write `is_unaligned`, `is_sorted` and `mixed_sq_headers`
# out so that they can be set as env.
# variables and handled downstream.
sys.stdout.write(
f"IS_UNALIGNED={int(is_unaligned)};" +
f"MIXED_HEADERS={int(mixed_headers)};" +
f"IS_SORTED={int(sorted_xam)}"
f"MIXED_SQ_HEADERS={int(mixed_sq_headers)};" +
f"IS_SORTED={int(sorted_xam)};" +
f"HAS_READS={int(has_reads)};"
)
logger.info(f"Checked (u)BAM headers in '{args.input_path}'.")
@ -68,4 +166,9 @@ def argparser():
"""Argument parser for entrypoint."""
parser = wf_parser("check_bam_headers_in_dir")
parser.add_argument("input_path", type=Path, help="Path to target directory")
parser.add_argument(
"--ref", type=Path, help="Optional reference FASTA file", required=False)
parser.add_argument(
"--ref_idx", type=Path,
help="Optional reference index FASTA file required for pysam", required=False)
return parser

48
bin/workflow_glue_r/R/bambu.R
View File

@ -1,8 +1,7 @@
bambu_arg_parser <- function() {
parser <- argparser::arg_parser("Run bambu transcript discovery and quantification.")
parser <- argparser::add_argument(parser, "--bam_dir", help = "Directory containing BAM files.")
parser <- argparser::add_argument(parser, "--bam_path", help = "Path to a single BAM file.")
parser <- argparser::add_argument(parser, "--sample_alias", help = "Alias to use for a single BAM file.")
parser <- argparser::add_argument(parser, "--bams", help = "Comma-separated BAM paths.")
parser <- argparser::add_argument(parser, "--aliases", help = "Comma-separated aliases for --bams.")
parser <- argparser::add_argument(parser, "--sample_sheet", help = "Optional sample sheet CSV.")
parser <- argparser::add_argument(parser, "--annotation", help = "Reference annotation GTF/GFF.")
parser <- argparser::add_argument(parser, "--genome", help = "Reference genome FASTA.")
@ -31,8 +30,11 @@ bambu_arg_parser <- function() {
bambu_validate_args <- function(argv) {
workflow_glue_r_require_args(argv, c("annotation", "genome", "out_dir"))
if (workflow_glue_r_arg_missing(argv$bam_dir) == workflow_glue_r_arg_missing(argv$bam_path)) {
stop("Provide exactly one of --bam_dir or --bam_path.", call. = FALSE)
if (workflow_glue_r_arg_missing(argv$bams)) {
stop("Missing required arguments: --bams", call. = FALSE)
}
if (workflow_glue_r_arg_missing(argv$aliases)) {
stop("Missing required arguments: --aliases", call. = FALSE)
}
if (!argv$transcriptome_mode %in% c("discover", "fixed_annotation")) {
@ -52,16 +54,9 @@ bambu_validate_args <- function(argv) {
invisible(argv)
}
bambu_strip_alias <- function(path) {
name <- basename(path)
name <- sub("\\.aligned\\.sorted\\.bam$", "", name)
tools::file_path_sans_ext(name)
}
bambu_resolve_inputs <- function(
argv,
bamfile_list_ctor = Rsamtools::BamFileList,
list_files_fn = base::list.files
bamfile_list_ctor = Rsamtools::BamFileList
) {
sample_df <- NULL
if (!workflow_glue_r_arg_missing(argv$sample_sheet)) {
@ -81,19 +76,14 @@ bambu_resolve_inputs <- function(
}
}
if (!workflow_glue_r_arg_missing(argv$bam_dir)) {
bam_paths <- sort(list_files_fn(argv$bam_dir, pattern = "\\.bam$", full.names = TRUE))
if (length(bam_paths) < 1) {
stop("No BAM files were found in bam_dir.", call. = FALSE)
}
aliases <- unname(vapply(bam_paths, bambu_strip_alias, character(1)))
} else {
bam_paths <- argv$bam_path
aliases <- if (!workflow_glue_r_arg_missing(argv$sample_alias)) {
argv$sample_alias
} else {
bambu_strip_alias(argv$bam_path)
}
bam_paths <- workflow_glue_r_parse_csv_list(argv$bams)
aliases <- workflow_glue_r_parse_csv_list(argv$aliases)
if (length(bam_paths) < 1) {
stop("No BAM files were provided in --bams.", call. = FALSE)
}
if (length(aliases) != length(bam_paths)) {
stop("Provide one alias per BAM in --bams.", call. = FALSE)
}
duplicate_bam_aliases <- unique(aliases[duplicated(aliases)])
@ -435,8 +425,7 @@ main_run_bambu <- function(
prepare_annotations_fn = bambu::prepareAnnotations,
gene_expression_fn = bambu::transcriptToGeneExpression,
write_gtf_fn = bambu::writeToGTF,
bamfile_list_ctor = Rsamtools::BamFileList,
list_files_fn = base::list.files
bamfile_list_ctor = Rsamtools::BamFileList
) {
set.seed(42)
suppressPackageStartupMessages({
@ -449,8 +438,7 @@ main_run_bambu <- function(
inputs <- bambu_resolve_inputs(
argv,
bamfile_list_ctor = bamfile_list_ctor,
list_files_fn = list_files_fn
bamfile_list_ctor = bamfile_list_ctor
)
annotation_obj <- prepare_annotations_fn(argv$annotation)
ndr_value <- bambu_resolve_ndr(argv)

76
bin/workflow_glue_r/R/de_analysis.R
View File

@ -149,20 +149,38 @@ de_extract_disp_gene_est <- function(object) {
S4Vectors::mcols(object)$dispGeneEst
}
# DESeq2's default geometric-mean size-factor estimator is undefined when
# every gene has at least one zero across samples, so use poscounts then.
de_choose_size_factor_type <- function(count_mat, context_label = "Count matrix") {
if (all(rowSums(count_mat == 0) > 0)) {
warning(
context_label,
" has every gene containing at least one zero; using DESeq2 size-factor estimation with sfType='poscounts'."
)
return("poscounts")
}
"ratio"
}
de_run_deseq_with_fallback <- function(
dds,
contrast_name,
out_dir
) {
sf_type <- de_choose_size_factor_type(
DESeq2::counts(dds),
context_label = paste0("DGE count matrix for ", contrast_name)
)
fallback_info <- list(
applied = FALSE,
method_used = "parametric",
reason = NULL,
diagnostic_file = NULL
diagnostic_file = NULL,
size_factor_type = sf_type
)
de_out <- tryCatch(
DESeq2::DESeq(dds, quiet = TRUE),
DESeq2::DESeq(dds, quiet = TRUE, sfType = sf_type),
error = function(err) {
if (!grepl(
"all gene-wise dispersion estimates are within 2 orders of magnitude",
@ -184,7 +202,7 @@ de_run_deseq_with_fallback <- function(
"Results will have reduced power and wider confidence intervals."
)
dds <- DESeq2::estimateSizeFactors(dds)
dds <- DESeq2::estimateSizeFactors(dds, type = sf_type)
dds <- DESeq2::estimateDispersionsGeneEst(dds)
dds <- de_set_dispersions(dds, de_extract_disp_gene_est(dds))
@ -242,7 +260,8 @@ de_run_deseq_with_fallback <- function(
applied = TRUE,
method_used = "gene-wise",
reason = conditionMessage(err),
diagnostic_file = basename(diag_file)
diagnostic_file = basename(diag_file),
size_factor_type = sf_type
)
DESeq2::nbinomWaldTest(dds)
@ -391,6 +410,17 @@ de_run_deseq2_result <- function(
)
}
de_estimate_size_factors_for_dexseq <- function(dxd, count_mat, coldata, sf_type) {
sf_dds <- DESeq2::DESeqDataSetFromMatrix(
countData = round(count_mat),
colData = coldata,
design = ~ 1
)
sf_dds <- DESeq2::estimateSizeFactors(sf_dds, type = sf_type)
DESeq2::sizeFactors(dxd) <- DESeq2::sizeFactors(sf_dds)
dxd
}
de_run_dexseq_result <- function(
tx_counts,
tx_meta,
@ -424,7 +454,16 @@ de_run_dexseq_result <- function(
featureID = tx_meta$TXNAME,
groupID = tx_meta$GENEID
)
dxd <- DESeq2::estimateSizeFactors(dxd)
dexseq_sf_type <- de_choose_size_factor_type(
round(tx_counts),
context_label = "DEXSeq transcript count matrix"
)
dxd <- de_estimate_size_factors_for_dexseq(
dxd,
tx_counts,
coldata,
dexseq_sf_type
)
dispersion_result <- de_estimate_dispersions_with_fallback(
dxd,
"DEXSeq",
@ -446,7 +485,8 @@ de_run_dexseq_result <- function(
dxd = dxd,
dxr = dxr,
dexseq_dispersion_method = dispersion_method,
dexseq_dispersion_reason = dispersion_reason
dexseq_dispersion_reason = dispersion_reason,
dexseq_size_factor_type = dexseq_sf_type
)
}, error = function(err) {
if (length(active_covariates) == 0 || !grepl(
@ -592,12 +632,14 @@ main_run_de_analysis <- function(argv) {
n_samples = nrow(contrast_samples),
n_target = sum(contrast_samples[[argv$condition_column]] == target_level),
n_reference = sum(contrast_samples[[argv$condition_column]] == reference_level),
deseq2_size_factor_method = "ratio",
deseq2_dispersion_fallback = list(
applied = FALSE,
method_used = "parametric",
reason = NULL,
diagnostic_file = NULL
),
dexseq_size_factor_method = "ratio",
dexseq_dispersion_method = "parametric",
dexseq_covariates_dropped = list()
)
@ -632,6 +674,9 @@ main_run_de_analysis <- function(argv) {
if (is.null(fallback_method) || identical(fallback_method, "")) {
fallback_method <- if (fallback_applied) "gene-wise" else "parametric"
}
if (!is.null(fallback$size_factor_type) && !identical(fallback$size_factor_type, "")) {
contrast_qc$deseq2_size_factor_method <- fallback$size_factor_type
}
contrast_qc$deseq2_dispersion_fallback <- list(
applied = fallback_applied,
method_used = fallback_method,
@ -737,6 +782,9 @@ main_run_de_analysis <- function(argv) {
contrast_qc$dtu_significant_transcripts <- 0
contrast_qc$dtu_significant_genes <- 0
} else {
if (!is.null(dex_res$dexseq_size_factor_type) && !identical(dex_res$dexseq_size_factor_type, "")) {
contrast_qc$dexseq_size_factor_method <- dex_res$dexseq_size_factor_type
}
if (!is.null(dex_res$dexseq_dispersion_method)) {
contrast_qc$dexseq_dispersion_method <- dex_res$dexseq_dispersion_method
}
@ -808,6 +856,7 @@ main_run_de_analysis <- function(argv) {
"",
"DGE Results:",
sprintf(" Genes tested: %d", contrast_qc$genes_tested),
sprintf(" Size factor method: %s", contrast_qc$deseq2_size_factor_method),
sprintf(" Significant (FDR < 0.05): %d", contrast_qc$dge_significant_fdr05),
sprintf(" Significant (FDR < 0.01): %d", contrast_qc$dge_significant_fdr01),
sprintf(" Upregulated: %d", contrast_qc$dge_upregulated),
@ -816,6 +865,7 @@ main_run_de_analysis <- function(argv) {
"DTU Results:",
sprintf(" Status: %s", contrast_qc$dtu_status),
sprintf(" Transcripts tested: %d", contrast_qc$transcripts_tested),
sprintf(" Size factor method: %s", contrast_qc$dexseq_size_factor_method),
if (contrast_qc$dtu_status == "SUCCESS") {
c(
sprintf(
@ -842,6 +892,10 @@ main_run_de_analysis <- function(argv) {
function(cqc) identical(cqc$deseq2_dispersion_fallback$method_used, "gene-wise"),
de_qc_stats$contrasts
))
deseq2_poscounts <- names(Filter(
function(cqc) identical(cqc$deseq2_size_factor_method, "poscounts"),
de_qc_stats$contrasts
))
dexseq_non_parametric <- names(Filter(
function(cqc) {
method <- cqc$dexseq_dispersion_method
@ -853,6 +907,10 @@ main_run_de_analysis <- function(argv) {
function(cqc) identical(cqc$dexseq_dispersion_method, "gene-wise"),
de_qc_stats$contrasts
))
dexseq_poscounts <- names(Filter(
function(cqc) identical(cqc$dexseq_size_factor_method, "poscounts"),
de_qc_stats$contrasts
))
dexseq_covariate_drop <- names(Filter(
function(cqc) length(cqc$dexseq_covariates_dropped) > 0,
de_qc_stats$contrasts
@ -867,10 +925,14 @@ main_run_de_analysis <- function(argv) {
deseq2_dispersion_fallback_contrast_names = as.list(deseq2_dispersion_fallbacks),
deseq2_gene_wise_contrasts = length(deseq2_gene_wise),
deseq2_gene_wise_contrast_names = as.list(deseq2_gene_wise),
deseq2_poscounts_contrasts = length(deseq2_poscounts),
deseq2_poscounts_contrast_names = as.list(deseq2_poscounts),
dexseq_non_parametric_dispersion_contrasts = length(dexseq_non_parametric),
dexseq_non_parametric_dispersion_contrast_names = as.list(dexseq_non_parametric),
dexseq_gene_wise_dispersion_contrasts = length(dexseq_gene_wise),
dexseq_gene_wise_dispersion_contrast_names = as.list(dexseq_gene_wise),
dexseq_poscounts_contrasts = length(dexseq_poscounts),
dexseq_poscounts_contrast_names = as.list(dexseq_poscounts),
dexseq_covariate_drop_contrasts = length(dexseq_covariate_drop),
dexseq_covariate_drop_contrast_names = as.list(dexseq_covariate_drop),
total_covariates_dropped = total_covariates_dropped
@ -910,7 +972,9 @@ main_run_de_analysis <- function(argv) {
sprintf(" %s:", cname),
sprintf(" Samples: %d (%d vs %d)", cqc$n_samples, cqc$n_target, cqc$n_reference),
sprintf(" DGE significant: %d genes (FDR<0.05)", cqc$dge_significant_fdr05),
sprintf(" DESeq2 size factors: %s", cqc$deseq2_size_factor_method),
sprintf(" DTU status: %s", cqc$dtu_status),
sprintf(" DEXSeq size factors: %s", cqc$dexseq_size_factor_method),
if (cqc$dtu_status == "SUCCESS") sprintf(" DTU significant: %d genes", cqc$dtu_significant_genes) else NULL
)
}),

109
bin/workflow_glue_r/tests/testthat/test_bambu.R
View File

@ -6,32 +6,31 @@
#' (bin/workflow_glue/prepare_annotation_reference.py) with pytest coverage.
#' These tests focus on bambu-specific validation and integration.
# Workflow must accept exactly one of --bam_dir or --bam_path, not both, not neither.
# Workflow requires explicit BAM paths and aliases.
# Fail fast with clear error rather than passing invalid inputs to bambu.
testthat::test_that("exactly one BAM input required", {
testthat::test_that("BAM inputs required", {
args <- list(
annotation = "annotation.gtf",
genome = "genome.fa",
out_dir = tempfile("bambu-out-"),
bam_dir = NULL,
bam_path = NULL,
bams = NULL,
aliases = NULL,
transcriptome_mode = "discover",
ndr = NULL
)
testthat::expect_error(
bambu_validate_args(args),
"Provide exactly one of --bam_dir or --bam_path"
"Missing required arguments: --bams"
)
args$bam_dir <- tempfile("bam-dir-")
args$bam_path <- "sample.bam"
args$bams <- "sampleA.bam"
testthat::expect_error(
bambu_validate_args(args),
"Provide exactly one of --bam_dir or --bam_path"
"Missing required arguments: --aliases"
)
args$bam_path <- NULL
args$aliases <- "sampleA"
testthat::expect_silent(bambu_validate_args(args))
})
@ -42,8 +41,8 @@ testthat::test_that("invalid discovery settings rejected", {
annotation = "annotation.gtf",
genome = "genome.fa",
out_dir = tempfile("bambu-out-"),
bam_dir = tempfile("bam-dir-"),
bam_path = NULL,
bams = "sampleA.bam",
aliases = "sampleA",
transcriptome_mode = "novel",
ndr = NULL
)
@ -72,48 +71,26 @@ testthat::test_that("invalid discovery settings rejected", {
testthat::expect_silent(bambu_validate_args(args))
})
# Derive clean sample aliases from BAM filenames by removing workflow-specific suffixes.
testthat::test_that("BAM suffixes stripped from aliases", {
testthat::expect_equal(
bambu_strip_alias("/tmp/sample.aligned.sorted.bam"),
"sample"
)
testthat::expect_equal(
bambu_strip_alias("/tmp/sample.bam"),
"sample"
)
})
# Fail fast if --bam_dir contains no BAM files rather than passing empty input to bambu.
testthat::test_that("empty BAM directory rejected", {
bam_dir <- tempfile("empty-bam-dir-")
dir.create(bam_dir)
# Fail fast if --bams is empty rather than passing empty input to bambu.
testthat::test_that("empty BAM list rejected", {
args <- list(
bam_dir = bam_dir,
bam_path = NULL,
sample_alias = NULL,
bams = "",
aliases = "",
sample_sheet = NULL
)
testthat::expect_error(
bambu_resolve_inputs(args, bamfile_list_ctor = function(paths, yieldSize) paths),
"No BAM files were found in bam_dir"
"No BAM files were provided in --bams"
)
})
# Sample aliases must be unique and explicitly specified.
# Derived aliases (from filenames) can collide: sample.bam vs sample.aligned.sorted.bam.
# Sample sheets must have an 'alias' column with unique values that match BAM inputs.
testthat::test_that("unique sample aliases required", {
bam_dir <- tempfile("ambiguous-bam-dir-")
dir.create(bam_dir)
file.create(file.path(bam_dir, "sample.bam"))
file.create(file.path(bam_dir, "sample.aligned.sorted.bam"))
args <- list(
bam_dir = bam_dir,
bam_path = NULL,
sample_alias = NULL,
bams = "sampleA.bam,sampleB.bam",
aliases = "sampleA,sampleA",
sample_sheet = NULL
)
testthat::expect_error(
@ -121,8 +98,11 @@ testthat::test_that("unique sample aliases required", {
"BAM aliases must be unique"
)
sample_bam <- tempfile(fileext = ".bam")
file.create(sample_bam)
args$aliases <- "sampleA"
testthat::expect_error(
bambu_resolve_inputs(args, bamfile_list_ctor = function(paths, yieldSize) paths),
"Provide one alias per BAM in --bams"
)
missing_alias_sheet <- tempfile(fileext = ".csv")
writeLines(
paste(
@ -133,9 +113,8 @@ testthat::test_that("unique sample aliases required", {
missing_alias_sheet
)
args <- list(
bam_dir = NULL,
bam_path = sample_bam,
sample_alias = "sampleA",
bams = "sampleA.bam",
aliases = "sampleA",
sample_sheet = missing_alias_sheet
)
testthat::expect_error(
@ -161,14 +140,9 @@ testthat::test_that("unique sample aliases required", {
})
# Sample sheet rows must align with BAM file order.
# If sheet lists samples in different order than filesystem, reorder sheet to match.
# If sheet lists samples in different order than CLI aliases, reorder sheet to match.
# If sheet is missing aliases found in BAMs, fail.
testthat::test_that("sample sheet reordered to match BAMs", {
bam_dir <- tempfile("bam-dir-")
dir.create(bam_dir)
file.create(file.path(bam_dir, "sampleA.aligned.sorted.bam"))
file.create(file.path(bam_dir, "sampleB.bam"))
sample_sheet <- tempfile(fileext = ".csv")
writeLines(
paste(
@ -181,9 +155,8 @@ testthat::test_that("sample sheet reordered to match BAMs", {
)
args <- list(
bam_dir = bam_dir,
bam_path = NULL,
sample_alias = NULL,
bams = "sampleA.aligned.sorted.bam,sampleB.bam",
aliases = "sampleA,sampleB",
sample_sheet = sample_sheet
)
resolved <- bambu_resolve_inputs(
@ -239,8 +212,8 @@ testthat::test_that("transcriptome mode mapped to bambu args", {
})
# End-to-end unit test with mocked bambu analysis function.
# Verifies --bam_dir input resolution, sample sheet reordering, and discovery settings.
testthat::test_that("bam_dir input with discovery mode", {
# Verifies --bams input resolution, sample sheet reordering, and discovery settings.
testthat::test_that("bams input with discovery mode", {
fixture_dir <- tempfile("bambu-discover-")
dir.create(fixture_dir)
bam_dir <- file.path(fixture_dir, "bams")
@ -281,9 +254,8 @@ testthat::test_that("bam_dir input with discovery mode", {
annotation = "annotation.gtf",
genome = "genome.fa",
out_dir = file.path(fixture_dir, "out"),
bam_dir = bam_dir,
bam_path = NULL,
sample_alias = NULL,
bams = paste(c(sample_a, sample_b), collapse = ","),
aliases = "sampleA,sampleB",
sample_sheet = sample_sheet,
transcriptome_mode = "discover",
ndr = 0.25,
@ -492,10 +464,10 @@ testthat::test_that("QC stats match filtered results", {
#
# End-to-end tests with real bambu library (not mocked):
# - Build BAMs from committed fixtures (reference.fa, annotation.gtf, reads.fastq)
# - Run `supeRglue bambu` with both --bam_path and --bam_dir inputs
# - Run `supeRglue bambu` with --bams input
# - Verify output files exist and contain data for downstream workflow steps
testthat::test_that("CLI single BAM input with fixed annotation", {
testthat::test_that("CLI single BAM in bams input with fixed annotation", {
fixture_dir <- tempfile("bambu-cli-")
dir.create(fixture_dir)
@ -510,8 +482,8 @@ testthat::test_that("CLI single BAM input with fixed annotation", {
"supeRglue",
c(
"bambu",
"--bam_path", bam_path,
"--sample_alias", "sampleA",
"--bams", bam_path,
"--aliases", "sampleA",
"--sample_sheet", sample_sheet,
"--annotation", annotation,
"--genome", reference,
@ -546,7 +518,7 @@ testthat::test_that("CLI single BAM input with fixed annotation", {
testthat::expect_gt(nrow(gene_counts), 0)
})
testthat::test_that("CLI bam_dir input preserves sample order", {
testthat::test_that("CLI bams input preserves sample order", {
fixture_dir <- tempfile("bambu-cli-dir-")
dir.create(fixture_dir)
bam_dir <- file.path(fixture_dir, "bams")
@ -573,7 +545,14 @@ testthat::test_that("CLI bam_dir input preserves sample order", {
"supeRglue",
c(
"bambu",
"--bam_dir", bam_dir,
"--bams", paste(
c(
file.path(bam_dir, "sampleA.aligned.sorted.bam"),
file.path(bam_dir, "sampleB.aligned.sorted.bam")
),
collapse = ","
),
"--aliases", "sampleA,sampleB",
"--sample_sheet", sample_sheet,
"--annotation", annotation,
"--genome", reference,

20
bin/workflow_glue_r/tests/testthat/test_de_analysis.R
View File

@ -388,6 +388,24 @@ testthat::test_that("underspecified designs rejected", {
# Fixture-driven tests for DESeq2/DEXSeq helper behaviour and metadata.
# These avoid dependency injection and exercise the real package code paths.
testthat::test_that("de_choose_size_factor_type uses poscounts for zero-heavy matrices", {
count_mat <- matrix(
c(
0, 4, 5,
3, 0, 2,
7, 1, 0
),
nrow = 3,
byrow = TRUE
)
testthat::expect_warning(
sf_type <- de_choose_size_factor_type(count_mat, "test matrix"),
"using DESeq2 size-factor estimation with sfType='poscounts'"
)
testthat::expect_equal(sf_type, "poscounts")
})
testthat::test_that("de_run_deseq_with_fallback returns structured metadata", {
testthat::skip_if_not_installed("DESeq2")
@ -721,7 +739,9 @@ testthat::test_that("CLI integration produces expected outputs", {
testthat::expect_true("analysis_fallbacks" %in% names(de_qc))
testthat::expect_true("contrasts" %in% names(de_qc))
contrast_qc <- de_qc$contrasts[["condition_treated_vs_control"]]
testthat::expect_true("deseq2_size_factor_method" %in% names(contrast_qc))
testthat::expect_true("deseq2_dispersion_fallback" %in% names(contrast_qc))
testthat::expect_true("dexseq_size_factor_method" %in% names(contrast_qc))
testthat::expect_true("dexseq_dispersion_method" %in% names(contrast_qc))
testthat::expect_true("dexseq_covariates_dropped" %in% names(contrast_qc))
})

1
docs/02_introduction.md
View File

@ -10,7 +10,6 @@ analysis.
The workflow supports:
+ transcript identification from either cDNA or direct RNA reads
+ optional cDNA preprocessing with `pychopper` before alignment
+ transcript discovery guided by a supplied genome and annotation
+ quantification against a supplied reference annotation
+ optional transcript classification and QC with `SQANTI3`

25
docs/06_pipeline_overview.md
View File

@ -60,15 +60,7 @@ barcode04,treated_rep2,test_sample,treated,b2
This example is suitable for a multiplexed run and also satisfies the minimum
requirements for a two-group DE/DTU comparison.
### 3. Optional cDNA preprocessing
When `--cdna_preprocess` is enabled for cDNA libraries, the workflow runs
`pychopper` before alignment to classify, orient, and trim full-length reads.
This preprocessing stage is controlled by `--cdna_kit`,
`--pychopper_backend`, and optional extra `--pychopper_opts`, and its outputs
are published alongside the ingress results for each sample.
### 4. Genome alignment
### 3. Genome alignment
Each sample is aligned to the supplied reference genome with
[`minimap2`](https://github.com/lh3/minimap2), then sorted and indexed with
@ -76,14 +68,14 @@ Each sample is aligned to the supplied reference genome with
`cohort/alignments/` are the main alignment files used for transcriptome
analysis, optional `SQANTI3` QC, and optional IGV viewing.
### 5. Cohort transcriptome construction
### 4. Cohort transcriptome construction
All aligned samples are analysed together with `bambu` to produce the primary
cohort transcriptome, transcript counts, gene counts, and the `RDS` objects used
for downstream differential analysis. This shared model is the main cohort-level
result and is published under `cohort/`.
### 6. Independent per-sample transcriptomes
### 5. Independent per-sample transcriptomes
Each sample is also processed separately with `bambu` so the workflow produces
sample-specific GTF, FASTA, count tables, and metadata under
@ -91,7 +83,7 @@ sample-specific GTF, FASTA, count tables, and metadata under
specific transcript models without changing the shared cohort transcriptome used
for DE/DTU.
### 7. Transcript sequence generation and QC
### 6. Transcript sequence generation and QC
Transcript FASTA files are derived from GTF plus genome using `gffread`.
When `--skip_sqanti` is not set, `SQANTI3` classifies the cohort and per-sample
@ -99,7 +91,7 @@ transcriptomes and produces structural QC summaries. The cohort `SQANTI3`
results live under `cohort/sqanti_cohort/`, while per-sample `SQANTI3`
directories are published under `samples/<alias>/<alias>_sqanti/`.
### 8. Optional DE and DTU analysis
### 7. Optional DE and DTU analysis
When `--de_analysis` is enabled, the workflow checks the experimental design,
runs `DESeq2` for differential gene expression, and runs `DEXSeq` for
@ -107,7 +99,7 @@ differential transcript usage. These analyses use the shared `bambu` outputs
and the design columns in the sample sheet, and each comparison is written to
its own subdirectory under `de_analysis/<contrast>/`.
### 9. What you need to provide
### 8. What you need to provide
The workflow's analysis is controlled by a user provided genome, annotation, and
`bambu` mode.
@ -119,13 +111,10 @@ The workflow's analysis is controlled by a user provided genome, annotation, and
* `--ref_transcriptome` has been removed; if you want annotation-based
quantification, use `--transcriptome_mode fixed_annotation` together with
`--ref_genome` and `--ref_annotation`
* `--cdna_preprocess` enables the `pychopper` cDNA preprocessing stage, and
`--cdna_kit`, `--pychopper_backend`, and `--pychopper_opts` control that stage
* when `--de_analysis` is enabled, the sample sheet must contain `alias`, the
primary condition column, and any requested columns named in `--covariates`
* `--cdna_preprocess` must not be combined with `--direct_rna`
### 10. How to read the output folder
### 9. How to read the output folder
The published outputs are organised around a small number of top-level
directories:

4
docs/07_input_parameters.md
View File

@ -17,7 +17,6 @@
| ref_annotation | string | Reference transcript annotation in GTF or GFF format. | Required in both discover and fixed_annotation modes. | |
| transcriptome_mode | string | How bambu should prepare the transcriptome model. | Use discover for reference-guided transcript discovery and quantification, or fixed_annotation for quantification only against the supplied annotation. | discover |
| direct_rna | boolean | Set this for direct RNA sequencing libraries. | | False |
| cdna_preprocess | boolean | Apply the optional cDNA preprocessing stage before alignment. | This is only relevant for cDNA libraries and must not be combined with `--direct_rna`. | False |
### Sample Options
@ -53,9 +52,6 @@
| threads | integer | Thread count to use for the core workflow processes. | | 4 |
| minimap2_opts | string | Extra command-line options to pass to minimap2. | | |
| ndr | number | Optional bambu novel discovery rate override. | | |
| cdna_kit | string | ONT cDNA kit identifier used for pychopper preprocessing. | Only used when `--cdna_preprocess` is enabled. The workflow derives the pychopper kit code from this value. | SQK-PCS109 |
| pychopper_backend | string | Primer-detection backend to use for pychopper preprocessing. | | edlib |
| pychopper_opts | string | Extra command-line options to pass to pychopper. | | |
| skip_sqanti | boolean | Skip SQANTI3 transcript classification and QC. | | False |
| sqanti_skip_orf | boolean | Skip ORF prediction during SQANTI3 QC. | | True |
| sqanti_extra_args | string | Extra command-line options to pass to SQANTI3. | | |

6
docs/08_outputs.md
View File

@ -7,9 +7,9 @@ Output files may be aggregated including information for all samples or provided
| Per-read stats | ingress_results/{{ alias }}/fastcat_stats/per-read-stats.tsv.gz | Read statistics for individual reads in a sample, when this output is enabled. | per-sample |
| Ingress reads | ingress_results/{{ alias }}/seqs.fastq.gz | Reads prepared from the input data for downstream analysis. | per-sample |
| Ingress metadata | ingress_results/{{ alias }}/metamap.json | Per-sample metadata used by the workflow. | per-sample |
| Aligned BAM | cohort/alignments/{{ alias }}.aligned.sorted.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
| Aligned BAM index | cohort/alignments/{{ alias }}.aligned.sorted.bam.bai | Index for the aligned BAM. | per-sample |
| Alignment summary | cohort/alignments/{{ alias }}.flagstat.txt | samtools flagstat output for the aligned BAM. | per-sample |
| Aligned BAM | cohort/alignments/{{ alias }}/reads.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
| Aligned BAM index | cohort/alignments/{{ alias }}/reads.bam.bai | Index for the aligned BAM. | per-sample |
| Alignment summary | cohort/alignments/{{ alias }}/bamstats.flagstat.tsv | bamstats flagstat summary for the aligned BAM. | per-sample |
| Reference and annotation preparation summary | cohort/reference/annotation_reference_summary.json | Summary of reference and annotation preparation, including seqname overlap, build/provider hints, and excluded unstranded annotation counts. | aggregated |
| Excluded unstranded annotation records | cohort/reference/unstranded_annotation.gtf | Full set of annotation records excluded because their strand was not '+' or '-'. Present only when unstranded records are found. | aggregated |
| Cohort transcriptome GTF | cohort/transcripts.gtf | Joint bambu transcript model used as the primary cohort transcriptome. | aggregated |

21
docs/10_troubleshooting.md
View File

@ -5,7 +5,6 @@
primary condition column, and any columns named in `--covariates`.
+ DE/DTU requires at least two condition levels and at least two samples per
level.
+ `--cdna_preprocess` must not be combined with `--direct_rna`.
+ See how to interpret common Nextflow exit codes
[here](https://labs.epi2me.io/trouble-shooting/).
@ -40,17 +39,6 @@ setting that chose how the workflow behaved. In the current version,
The options `--transcriptome_mode discover` or `--transcriptome_mode
fixed_annotation` should be used to choose between the modes of operation.
#### I expected pychopper-style preprocessing for cDNA libraries
cDNA preprocessing is available through
`--cdna_preprocess`, which runs `pychopper` before alignment. The
parameters for that stage are `--cdna_kit`, `--pychopper_backend`, and
`--pychopper_opts`.
Enable `--cdna_preprocess` for cDNA libraries when you want
`pychopper` preprocessing, or use `--direct_rna` without
`--cdna_preprocess` for direct RNA data.
#### I cannot find the old flat DE output files
In the previous workflow version, DE files appeared directly under `de_analysis/`.
@ -67,9 +55,7 @@ In the current versions, the output folder is organised around `ingress_results/
`igv_reference/`.
Primary shared transcriptome results are under `cohort/`, and
`samples/<alias>/` for sample-specific models. If `--cdna_preprocess` is enabled,
look under `ingress_results/<alias>/` for the `pychopper` outputs associated
with each sample.
`samples/<alias>/` for sample-specific models.
#### My DE/DTU run fails because of `--sample_sheet`
@ -83,11 +69,6 @@ are required, and each level must contain at least two samples.
What to change: verify the sample sheet columns first, then check
`--condition_column`, `--covariates`, and `--reference_level`.
#### I combined `--direct_rna` with `--cdna_preprocess`
This is not a valid combination of parameters. Use `--direct_rna` for direct RNA
libraries and leave `--cdna_preprocess` unset.
#### I hit genome/annotation validation or strand-related annotation warnings
The workflow validates that the annotation and genome share

76
lib/common.nf
View File

@ -61,3 +61,79 @@ process configure_igv {
"""
}
process minimap2_alignment {
label "wf_common"
cpus { bsargs["alignment_threads"] }
memory { bsargs["minimap2_memory"][task.attempt - 1] }
maxRetries { bsargs["minimap2_memory"].size() - 1 }
errorStrategy = {task.exitStatus in [137,140] ? 'retry' : 'finish'}
input:
tuple val(meta), path(reads, stageAs: 'reads/*'), path(reference), path(ref_index)
tuple val(align_ext), val(index_ext) // either [bam, bai] or [cram, crai]
val bsargs
output:
tuple val(meta),
path("reads.${align_ext}"),
path("reads.${align_ext}.${index_ext}"),
path("bamstats_results"),
emit: alignment
script:
String reset_cmd_body = "samtools reset -x tp,cm,s1,s2,NM,MD,AS,SA,ms,nn,ts,cg,cs,dv,de,rl"
String fastq_cmd_body = "samtools fastq -T *"
// Default required threads is 6
// Samtools x3 and bamstats will all be single threaded
def minimap2_threads = Math.max(task.cpus - 4 , 1)
def per_read_stats_arg = bsargs["per_read_stats"] ? "| bgzip > bamstats_results/bamstats.readstats.tsv.gz" : " > /dev/null"
def bam_input_cmd = """samtools view -H --no-PG reads/\"\$(ls reads | head -n 1)\" \
> reads.header && samtools cat reads/* \
| ${reset_cmd_body} --no-PG - -o - \
| ${fastq_cmd_body} - """
def fastq_input_cmd = "touch reads.header && fastcat --reheader reads/*"
def bam_or_fastq_input = params.bam ? bam_input_cmd : fastq_input_cmd
def minimap2_opts = bsargs["minimap2_opts"]
"""
rm -rf bamstats_results
mkdir bamstats_results
${bam_or_fastq_input} \
| minimap2 -y -t ${minimap2_threads} -a ${minimap2_opts} --cap-kalloc 100m --cap-sw-mem 50m \
${reference} - \
| workflow-glue reheader_samstream reads.header \
--insert \$'@PG\\tID:reset\\tPN:samtools\\tCL:${reset_cmd_body}' \
--insert \$'@PG\\tID:fastq\\tPN:samtools\\tCL:${fastq_cmd_body}' \
| samtools sort -u -O BAM - \
| tee >(samtools view --reference ${reference} \
--write-index -o reads.${align_ext}##idx##reads.${align_ext}.${index_ext} -) \
| bamstats -s ${meta.alias} -u \
-f bamstats_results/bamstats.flagstat.tsv \
-i bamstats_results/bamstats.runids.tsv \
-l bamstats_results/bamstats.basecallers.tsv \
--histograms histograms - \
${per_read_stats_arg}
# post-pipe bamstats tidying
mv histograms/* bamstats_results/
# get n_seqs from flagstats - need to sum them up
awk 'NR==1{for (i=1; i<=NF; i++) {ix[\$i] = i}} NR>1 {c+=\$ix["total"]} END{print c}' \
bamstats_results/bamstats.flagstat.tsv > bamstats_results/n_seqs
# get unique run IDs (we add `-F '\\t'` as `awk` uses any stretch of whitespace
# as field delimiter otherwise and thus ignore empty columns)
awk -F '\\t' '
NR==1 {for (i=1; i<=NF; i++) {ix[\$i] = i}}
# only print run_id if present
NR>1 && \$ix["run_id"] != "" {print \$ix["run_id"]}
' bamstats_results/bamstats.runids.tsv | sort | uniq > bamstats_results/run_ids
# get unique basecall models
awk -F '\\t' '
NR==1 {for (i=1; i<=NF; i++) {ix[\$i] = i}}
# only print model if present
NR>1 && \$ix["basecaller"] != "" {print \$ix["basecaller"]}
' bamstats_results/bamstats.basecallers.tsv | sort | uniq > bamstats_results/basecallers ;
"""
}

239
lib/ingress.nf
View File

@ -4,6 +4,10 @@ import ArgumentParser
N_OPEN_FILES_LIMIT = 128
OPTIONAL_FILE = file("$projectDir/data/OPTIONAL_FILE")
include { minimap2_alignment as bam_alignment;
minimap2_alignment as fastq_alignment; } from './common'
/**
* Check if a file ends with one of the target extensions.
@ -171,6 +175,14 @@ def add_number_of_reads_to_meta(ch, String input_type_format) {
* - "allow_multiple_basecall_models": emit data of samples that had more than one
* basecall model; if this is `false`, such samples will be emitted as `[meta, null,
* null]`
* - "minimap2_memory": A list of memory options to pass to minimap2. Alignment will be retried
* with next in list if it fails due to memory issues. (default: ["8G", "15G", "31G"])
* - "minimap2_opts": string with options to pass to minimap2 (default: "-x lr:hq")
* - "alignment_threads": number of threads to use for alignment process (default: 6)
* - "output_xam_fmt": alignment output format, `bam` outputs a BAM file with index (.bam, .bai),
`cram` outputs a CRAM file with index (.cram, .crai)
* @param aln_ref_ch: optional channel with a reference tuple (ref, ref_idx) to align against.
* If provided alignment will be attempted if inputs are unaligned or not already aligned to this ref.
* @return: channel of `[Map(alias, barcode, type, ...), Path|null, Path|null]`.
* The first element is a map with metadata, the second is the path to the
* `.fastq.gz` file with the (potentially concatenated) sequences and the third is
@ -178,7 +190,7 @@ def add_number_of_reads_to_meta(ch, String input_type_format) {
* `null` for sample sheet entries for which no corresponding barcode directory was
* found. The third element is `null` if there were no reads.
*/
def fastq_ingress(Map arguments)
def fastq_ingress(Map arguments, aln_ref_ch = null)
{
// check arguments
Map margs = parse_arguments(
@ -194,42 +206,99 @@ def fastq_ingress(Map arguments)
def input = get_valid_inputs(margs, fq_extensions)
def ch_result = fastcat(input.files.mix(input.dirs), margs, "FASTQ")
// TODO: xam_ingress mixes in a .no_files channel here. Do we need to do the same?
// Alignment with reference if provided
if (aln_ref_ch) {
// The above may have returned a channel with multiple fastqs if chunking
// is enabled. Flatten this and add a groupKey to meta information which
// states the number of sibling files. This can be later used as the key
// for .groupTuple() on a channel in order to get all results for a sample
// We don't decorate "alias" with a count because that messes up downstream
// serialisation.
// Mix in the missing files from the sample sheet
// Add in a unique key for every emission
def ch_spread_result = ch_result
.mix (input.missing.map { meta, files -> [meta, files, null] })
.map { meta, files, stats ->
// new `arity: '1..*'` would be nice here
files = files instanceof List ? files : [files]
def new_keys = [
"group_key": groupKey(meta["alias"], files.size()),
"n_fastq": files.size()]
def grp_index = (0..<files.size()).collect()
[meta + new_keys, files, grp_index, stats]
if (margs["output_xam_fmt"] == "bam"){
output_xam_fmt = ["bam", "bai"]
}
.transpose(by: [1, 2]) // spread multiple fastq files into separate emissions
.map { meta, files, grp_i, stats ->
def new_keys = [
"group_index": "${meta["alias"]}_${grp_i}"]
[meta + new_keys, files, stats]
else {
output_xam_fmt = ["cram", "crai"]
}
mm2_aln = fastq_alignment(
input.files.mix(input.dirs).combine(aln_ref_ch),
output_xam_fmt,
margs
)
// Add back samples from sample sheet with missing barcode dir
aligned_and_missing = mm2_aln.mix(
input.missing.map { meta, files -> [meta, files, null, null] }
)
// add number of reads to meta
updated_ch = add_number_of_reads_to_meta(
aligned_and_missing
| map { meta, bam, bai, stats ->
[meta, [bam, bai], stats]
},
"xam"
)
// add run IDs, and basecall models to meta
ch_result = add_run_IDs_and_basecall_models_to_meta(
updated_ch,
margs.allow_multiple_basecall_models
) | map { it ->
def (meta, bam, bai, stats) = it.flatten()
def xam = meta.src_xam
def xai = meta.src_xai
// S3 paths are remote references; only keep local file paths for downstream processing
if (meta.src_xam){
xam = meta.src_xam.startsWith('s3://') ? null : meta.src_xam
xai = meta.src_xam.startsWith('s3://') ? null : meta.src_xai
}
[ meta + [src_xam: xam, src_xai: xai], bam, bai, stats ]
}
// add number of reads, run IDs, and basecall models to meta
def ch_final = add_number_of_reads_to_meta(ch_spread_result, "fastq")
ch_final = add_run_IDs_and_basecall_models_to_meta(
ch_final, margs.allow_multiple_basecall_models
)
return ch_final
return ch_result
} else {
// If no alignment ref is provided, just process FASTQ as usual.
def ch_fastcat = fastcat(input.files.mix(input.dirs), margs, "FASTQ")
def ch_fastcat_branched = ch_fastcat.branch { meta, files, stats ->
def file_list = files ? (files instanceof List ? files : [files]) : []
no_files: file_list.size() == 0
has_files: true
}
ch_result = ch_fastcat_branched.has_files
def ch_no_files = ch_fastcat_branched.no_files
.map { meta, files, stats -> [meta, null, null] }
// fastcat may have returned a channel with multiple fastqs if chunking
// is enabled. Flatten this and add a groupKey to meta information which
// states the number of sibling files. This can be later used as the key
// for .groupTuple() on a channel in order to get all results for a sample
// We don't decorate "alias" with a count because that messes up downstream
// serialisation.
// Mix in the missing files from the sample sheet
// Add in a unique key for every emission
ch_spread_result = ch_result
.mix(input.missing.map { meta, files -> [meta, files, null] })
.mix(ch_no_files)
.map { meta, files, stats ->
// new `arity: '1..*'` would be nice here
files = files instanceof List ? files : [files]
def new_keys = [
"group_key": groupKey(meta["alias"], files.size()),
"n_fastq": files.size()]
def grp_index = (0..<files.size()).collect()
[meta + new_keys, files, grp_index, stats]
}
.transpose(by: [1, 2]) // spread multiple fastq files into separate emissions
.map { meta, files, grp_i, stats ->
def new_keys = [
"group_index": "${meta["alias"]}_${grp_i}"]
[meta + new_keys, files, stats]
}
// add number of reads, run IDs, and basecall models to meta
ch_final = add_number_of_reads_to_meta(ch_spread_result, "fastq")
ch_final = add_run_IDs_and_basecall_models_to_meta(
ch_final, margs.allow_multiple_basecall_models
)
return ch_final
}
}
@ -251,7 +320,7 @@ def fastq_ingress(Map arguments)
* directories.
* - "keep_unaligned": boolean whether to include uBAM files
* - "return_fastq": boolean whether to convert to FASTQ (this will always run
* `fastcat`)
* `fastcat`). Cannot be set to true if an alignment reference is provided.
* - "fastcat_extra_args": string with extra arguments to pass to `fastcat`
* - "required_sample_types": list of zero or more required sample types expected to be present
* in the sample sheet
@ -260,7 +329,15 @@ def fastq_ingress(Map arguments)
* - "fastq_chunk": null or a number of reads to place into chunked FASTQ files
* - "allow_multiple_basecall_models": boolean. If true, emit data of samples that had more than one
* basecall model; if this is `false`, such samples will be emitted as `[meta, null,
* null]`
null]`
* - "minimap2_memory": a list of memory options to pass to minimap2, alignment will be retried
* with next in list if it fails due to memory issues. (default: ["8G", "15G", "31G"])
* - "minimap2_opts": string with options to pass to minimap2 (default: "-x lr:hq")
* - "alignment_threads": number of threads to use for alignment process (default: 6)
* - "output_xam_fmt": alignment output format, `bam` outputs an BAM file with index (.bam, .bai),
`cram` outputs a CRAM file with index (.cram, .crai)
* @param aln_ref_ch: optional channel with a reference tuple (ref, ref_idx) to align against.
* If provided alignment will be attempted if inputs are unaligned or not already aligned to this ref.
* @return: channel of `[Map(alias, barcode, type, ...), Path|null, Path|null]`.
* The first element is a map with metadata, the second is the path to the
* `.bam` file with the (potentially merged) sequences and the third is
@ -269,7 +346,7 @@ def fastq_ingress(Map arguments)
* found and for samples with only uBAM files when `keep_unaligned: false`. The third
* element is `null` if `bamstats` was not run.
*/
def xam_ingress(Map arguments)
def xam_ingress(Map arguments, aln_ref_ch = null)
{
// check arguments
Map margs = parse_arguments(
@ -283,13 +360,25 @@ def xam_ingress(Map arguments)
)
margs["fastq_chunk"] ?= 0 // cant pass null through channel
if (margs["return_fastq"] && aln_ref_ch){
error "`return_fastq` ingress argument cannot be true when alignment is enabled. Alignment always produces XAM-format outputs."
}
// we only accept BAM or uBAM for now (i.e. no SAM or CRAM)
ArrayList xam_extensions = [".bam", ".ubam"]
def input = get_valid_inputs(margs, xam_extensions)
if (aln_ref_ch) {
aln_ref = aln_ref_ch
} else {
aln_ref = Channel.of(
tuple(OPTIONAL_FILE, OPTIONAL_FILE)
)
}
// check BAM headers to see if any samples are uBAM
ch_result = input.dirs
ch_check_bams = input.dirs
| map { meta, path -> [meta, get_target_files_in_dir(path, xam_extensions, margs)] }
| mix(input.files)
| map{
@ -313,21 +402,65 @@ def xam_ingress(Map arguments)
}
[meta + [src_xam: src_xam, src_xai: src_xai], paths]
}
| combine(aln_ref)
| checkBamHeaders
| map { meta, paths, is_unaligned_env, mixed_headers_env, is_sorted_env ->
| map { meta, paths, is_unaligned_env, mixed_sq_headers_env, is_sorted_env, has_reads_env ->
// convert the env. variables from strings ('0' or '1') into bools
boolean is_unaligned = is_unaligned_env as int as boolean
boolean mixed_headers = mixed_headers_env as int as boolean
boolean mixed_sq_headers = mixed_sq_headers_env as int as boolean
boolean is_sorted = is_sorted_env as int as boolean
// if no reads, no error (multisample may have some empty bams) but do not attempt alignment
boolean has_reads = has_reads_env as int as boolean
// throw an error if there was a sample with mixed headers
if (mixed_headers) {
if (mixed_sq_headers) {
error "Found mixed headers in (u)BAM files of sample '${meta.alias}'."
}
// add `is_unaligned` to the metamap (note the use of `+` to create a copy of
// `meta` to avoid modifying every item in the channel;
// https://github.com/nextflow-io/nextflow/issues/2660)
[meta + [is_unaligned: is_unaligned, is_sorted: is_sorted], paths]
[meta + [is_unaligned: is_unaligned, is_sorted: is_sorted, has_reads: has_reads], paths]
}
// Handle alignment
if (aln_ref_ch) {
alignment_fork = ch_check_bams
| branch {
meta, paths ->
to_align: (meta.is_unaligned == true) && (meta.has_reads == true)
noalign: true
}
if (margs["output_xam_fmt"] == "bam"){
output_xam_fmt = ["bam", "bai"]
}
else {
output_xam_fmt = ["cram", "crai"]
}
mm2_aln = bam_alignment(
alignment_fork.to_align.combine(aln_ref_ch),
output_xam_fmt,
margs
)
// Update meta is unaligned
mm2_aln_final = mm2_aln.alignment.map{
meta, xam, xai, stats ->
// remove has reads from meta as no longer required
def newmeta = meta.findAll { k, v -> k != 'has_reads' }
[newmeta + [is_unaligned: false], xam, xai, stats]
}
// Process BAM files that do not require realignment by passing them through the standard downstream steps (merging, sorting, indexing, etc.)
ch_result_tmp = alignment_fork.noalign.map{
meta, paths ->
def newmeta = meta.findAll { k, v -> k != 'has_reads' }
[newmeta, paths]
}
} else {
// If no alignment reference provided process all BAM's as usual.
ch_result_tmp = ch_check_bams.map{ meta, paths ->
def newmeta = meta.findAll { k, v -> k != 'has_reads' }
[newmeta, paths]}
}
ch_result = ch_result_tmp
| branch { meta, paths ->
// set `paths` to `null` for uBAM samples if unallowed (they will be added to
// the results channel in shape of `[meta, null]` at the end of the function
@ -485,7 +618,8 @@ def xam_ingress(Map arguments)
}
// Combine all possible inputs
ch_result = ch_missing | mix(
ch_result = ch_missing
| mix(
ch_no_op,
ch_indexed,
ch_merged,
@ -514,6 +648,11 @@ def xam_ingress(Map arguments)
ch_result.is_null.map{it + [null]}
)
// Add back aligned if alignment reference provided
if (aln_ref_ch) {
ch_result = ch_result.mix(mm2_aln_final)
}
// Remove metadata that are unnecessary downstream:
// meta.src_xai: not needed, as it will be part of the channel as a file
// meta.is_sorted: if data are aligned, they will also be sorted/indexed
@ -550,13 +689,14 @@ def xam_ingress(Map arguments)
meta, bam, bai, stats ->
def xam = meta.src_xam
def xai = meta.src_xai
// S3 paths are remote references;
// only keep local file paths for downstream processing
if (meta.src_xam){
xam = meta.src_xam.startsWith('s3://') ? null : meta.src_xam
xai = meta.src_xam.startsWith('s3://') ? null : meta.src_xai
}
[ meta + [src_xam: xam, src_xai: xai], bam, bai, stats ]
}
return ch_result
}
@ -629,18 +769,21 @@ process checkBamHeaders {
label "wf_common"
cpus 1
memory "2 GB"
input: tuple val(meta), path("input_dir/reads*.bam")
input: tuple val(meta), path("input_dir/reads*.bam"), path(ref, stageAs: "ref/*"),
path(ref_index, stageAs: "ref_index/*")
output:
tuple(
val(meta),
path("input_dir/reads*.bam", includeInputs: true),
env(IS_UNALIGNED),
env(MIXED_HEADERS),
env(MIXED_SQ_HEADERS),
env(IS_SORTED),
env(HAS_READS),
)
script:
String ref_arg = ref.fileName.name == OPTIONAL_FILE.name ? "" : "--ref $ref --ref_idx $ref_index"
"""
workflow-glue check_bam_headers_in_dir input_dir > env.vars
workflow-glue check_bam_headers_in_dir input_dir $ref_arg > env.vars
source env.vars
"""
}
@ -836,6 +979,10 @@ Map parse_arguments(String func_name, Map arguments, Map extra_kwargs=[:]) {
"required_sample_types": [],
"per_read_stats": false,
"allow_multiple_basecall_models": false,
"minimap2_memory": ["8GB", "15GB", "31GB"],
"minimap2_opts": "-x lr:hq",
"alignment_threads": 6,
"output_xam_fmt": "bam" // or cram
]
ArgumentParser parser = new ArgumentParser(
args: required_args,

246
lib/reference.nf Normal file
View File

@ -0,0 +1,246 @@
// This subworkflow pre-processes an input reference
// genome by preparing the appropriate indexes and cache.
// By default, the workflow takes an input reference,
// decompresses it if it is gzipped, and either imports
// an existing index (if available) or generates one
// using `samtools faidx`.
// If requested by the user, the workflow can also:
// 1. Generate a CRAM cache, using the `REF_PATH`
// environment variable. IMPORTANT: If the ingress workflow
// outputs CRAM files (when output_xam_fmt='cram'), a CRAM cache
// MUST be generated here and passed to downstream processes.
// Even though users cannot currently input CRAMs, the ingress
// workflow can generate them, and downstream processes require
// the cache to decompress CRAM files.
// 2. Generate a minimap2 `.mmi` index for faster alignment.
Map colors = NfcoreTemplate.logColours(params.monochrome_logs)
// Argument parser
Map parse_reference(Map arguments) {
ArgumentParser parser = new ArgumentParser(
args:[
],
kwargs:[
"output_cache": false,
"output_mmi": false,
"mmi_opts": "-x lr:hq",
"mmi_memory": ["8GB", "15GB", "31GB"],
],
name: "reference_ingress")
return parser.parse_args(arguments)
}
// Process to generate the CRAM cache and
// create the REF_PATH variable.
// NOTE: The CRAM cache is essential for any downstream processes that
// need to decompress or work with CRAM files. Even though users cannot
// currently input CRAMs, if output_xam_fmt='cram' is used in ingress,
// downstream processes will require this cache to properly handle the
// CRAM output from ingress.
process cram_cache {
label "wf_common"
cpus 1
memory 4.GB
input:
path reference
output:
tuple path("ref_cache/"), env(REF_PATH), emit: ref_cache
script:
"""
# Invoke from binary installed to container PATH
seq_cache_populate.pl -root ref_cache/ "${reference}"
REF_PATH="ref_cache/%2s/%2s/%s"
"""
}
// Process to create the faidx index
// main.nf currently handles publishing the faidx if required
// maintain naming so it matches input ref when output
process faidx {
label "wf_common"
cpus 1
memory 4.GB
input:
path(ref)
output:
path("${ref}.fai")
script:
"""
samtools faidx "${ref}"
"""
}
// Process to create the faidx indexes for a gzipped reference
process gz_faidx {
publishDir "${params.out_dir}", mode: 'copy', pattern: "*"
label "wf_common"
cpus 1
memory 4.GB
// If a user provides a non-bgzipped file indexing will fail, which should be tolerated.
// With error_strategy set to 'ignore', this will result in an empty channel
// that can be checked for downstream using `ifEmpty()`
errorStrategy 'ignore'
input:
path(ref)
output:
tuple path(ref), path("${ref}.fai"), path("${ref}.gzi")
script:
"""
samtools faidx "${ref}"
"""
}
// Decompress the reference genome
// NOTE -f required to decompress symlink
process decompress_ref {
label "wf_common"
cpus 1
memory 4.GB
input:
path "ref.fa.gz"
output:
path "ref.fa", emit: decompressed_ref
"""
gzip -df ref.fa.gz
"""
}
// Prepare minimap2 .mmi index
process make_mmi {
memory { margs["mmi_memory"][task.attempt - 1] }
maxRetries { margs["mmi_memory"].size() - 1 }
errorStrategy 'retry'
label "wf_common"
cpus 3
input:
path("ref.fa")
val margs
output:
path("ref.mmi")
script:
"""
minimap2 -t ${task.cpus} ${margs["mmi_opts"]} -d ref.mmi ref.fa
"""
}
// Workflow to prepare the reference genome and its indexes.
/**
* This workflow accepts either a string file path or a Nextflow channel as the reference input,
* making it flexible for both direct file inputs and outputs from upstream processes.
*
* IMPORTANT: If your workflow uses CRAM output from ingress (output_xam_fmt='cram'),
* you must set "output_cache": true to generate the CRAM cache. This cache is required
* by downstream processes to properly handle CRAM files, even though users cannot currently
* input CRAMs directly.
*
* inputs:
* - reference_path_or_channel: (String or Channel) Path to reference FASTA file or a Nextflow channel emitting a reference file
* - arguments: (Map) Optional configuration parameters:
* - "output_cache": (boolean) If true, generate a CRAM cache for the reference (default: false)
* REQUIRED if using CRAM output from ingress for downstream processing.
* - "output_mmi": (boolean) If true, generate a minimap2 .mmi index for the reference (default: false)
* - "mmi_opts": (string) Options string passed to minimap2 for index creation (default: "-x lr:hq")
* - "mmi_memory": (list) Memory options to pass to minimap2 index; creation will be retried with next in list if it fails due to memory issues (default: ["8GB", "15GB", "31GB"])
*
* emits:
* - "ref_tuple": Tuple of (reference, fai_index)
* - "ref_cache": CRAM cache for the reference genome if "output_cache" is true.
* Required for downstream processes when CRAM files are used.
* - "ref_mmi": Minimap2 index for the reference genome if "output_mmi" is true
* - "ref_gzidx": Gzip index for the reference genome (useful for IGV)
*/
workflow prepare_reference {
take:
reference_path_or_channel
arguments
main:
Map margs = parse_reference(arguments)
// Determine if input is a string path or channel
Boolean is_string = reference_path_or_channel in String
// Base ref channel
if (!is_string) {
// It's a channel - use it directly
ref = reference_path_or_channel
is_compressed = false // Assume uncompressed when coming from channel
gzindexes = Channel.empty()
} else {
// It's a string path
if (file(reference_path_or_channel).exists()){
ref = Channel.fromPath(reference_path_or_channel)
} else {
throw new Exception(colors.red + "File ${reference_path_or_channel} not found." + colors.reset)
}
// Check if compressed
is_compressed = reference_path_or_channel.toLowerCase().endsWith("gz")
if (is_compressed) {
// Define indexes names.
String input_fai_index = "${reference_path_or_channel}.fai"
String input_gzi_index = "${reference_path_or_channel}.gzi"
// Decompress reference genome.
ref = decompress_ref(ref)
// Check whether the input gzref is indexed. If so, pass these as indexes.
// Otherwise, generate the gzip + fai indexes for the compressed reference.
if (file(input_fai_index).exists() && file(input_gzi_index).exists()){
gzindexes = Channel.fromPath(reference_path_or_channel)
| mix(
Channel.fromPath(input_fai_index),
Channel.fromPath(input_gzi_index)
)
} else {
gzindexes = gz_faidx(Channel.fromPath(reference_path_or_channel))
gzindexes.ifEmpty{
if (params.containsKey("igv") && params.igv){
log.warn """\
The input reference is compressed but not with bgzip, which is required to create an index.
The workflow will proceed but it will not be possible to load the reference in the IGV Viewer.
To use the IGV Viewer, provide an uncompressed, or bgzip compressed version of the input reference next time you run the workflow.
""".stripIndent()
}
}
}
} else {
gzindexes = Channel.empty()
}
}
// Generate fai index if the file is either compressed, or if fai doesn't exists
if (!is_compressed) {
if (is_string && file("${reference_path_or_channel}.fai").exists()){
ref_idx = Channel.fromPath("${reference_path_or_channel}.fai")
} else {
ref_idx = faidx(ref)
}
} else {
ref_idx = faidx(ref)
}
// Combine ref and its fai index
ref_tuple = ref.combine(ref_idx)
// Generate CRAM cache
if (margs.output_cache){
cram_cache(ref)
ref_cache = cram_cache.out.ref_cache
} else {
ref_cache = null
}
// Generate mmi index
if (margs.output_mmi){
ref_mmi = make_mmi(ref, margs)
} else {
ref_mmi = null
}
emit:
ref_tuple = ref_tuple
ref_cache = ref_cache
ref_mmi = ref_mmi
ref_gzidx = gzindexes
}

254
main.nf
View File

@ -5,6 +5,7 @@ nextflow.enable.dsl = 2
include { fastq_ingress; xam_ingress } from './lib/ingress'
include { getParams; configure_igv } from './lib/common'
include { prepare_reference } from './lib/reference'
include { transcriptome_analysis } from './subworkflows/transcriptome'
include { differential_expression } from './subworkflows/differential_expression'
@ -29,69 +30,6 @@ process getVersions {
}
process collectIngressResultsInDir {
label "wf_common"
cpus 1
memory "2 GB"
input:
tuple val(meta), path(reads, stageAs: "reads/*"), path(stats, stageAs: "stats/*")
output:
path "out/*", emit: out
script:
String outdir = "out/${meta["alias"].replaceAll("'", "'\\\\''")}"
String meta_json = new JsonBuilder(meta).toPrettyString().replaceAll("'", "'\\\\''")
String stats_arg = stats.fileName.name == OPTIONAL_FILE.name ? "" : stats
def read_args = reads instanceof java.util.ArrayList ? reads.join(" ") : reads
"""
mkdir -p '${outdir}'
echo '${meta_json}' > metamap.json
mv metamap.json ${read_args} ${stats_arg} '${outdir}'
"""
}
process preprocess_reads {
label "wf_transcriptomes_pychopper"
cpus { params.threads ?: 4 }
memory "8 GB"
input:
tuple val(meta), path(reads, stageAs: "reads/*")
output:
tuple val(meta.alias), path("${meta.alias}_pychopper_output/${meta.alias}_full_length_reads.fastq"), emit: full_len_reads
tuple val(meta.alias), path("${meta.alias}_pychopper_output"), emit: dir
script:
String backend = params.pychopper_backend ?: "edlib"
String extra = params.pychopper_opts ?: ""
String cdna_kit = params.cdna_kit ? params.cdna_kit.tokenize("-")[-1] : ""
def read_args = reads instanceof java.util.ArrayList ? reads.join(" ") : reads
"""
cat ${read_args} > seqs.fastq.gz
pychopper -t ${task.cpus} -k ${cdna_kit} -m ${backend} ${extra} \
seqs.fastq.gz "${meta.alias}_full_length_reads.fastq"
awk '
BEGIN { FS = OFS = "\\t" }
NR == 1 {
for (i = 1; i <= NF; i++) {
idx[\$i] = i
}
print "Classification", "Value"
next
}
\$idx["Category"] == "Classification" {
print \$idx["Name"], \$idx["Value"]
}
' pychopper.tsv > pychopper_summary.tsv
mkdir -p "${meta.alias}_pychopper_output"
find . -maxdepth 1 -mindepth 1 \
! -name "seqs.fastq.gz" \
! -name "${meta.alias}_pychopper_output" \
-exec mv -t "${meta.alias}_pychopper_output" {} +
"""
}
process makeReport {
label "wf_common"
publishDir "${params.out_dir}", mode: 'copy', pattern: "wf-transcriptomes-report.html"
@ -101,10 +39,8 @@ process makeReport {
tuple val(metadata), path(stats, stageAs: "stats_*")
path "versions/*"
path "params.json"
path alignment_stats, stageAs: "alignment_stats/*"
path cohort_dir, stageAs: "cohort"
path sample_dirs, stageAs: "samples/*"
path pychopper_dirs, stageAs: "pychopper/*"
path sqanti_dirs, stageAs: "sqanti/*"
path de_files
val wf_version
@ -112,20 +48,17 @@ process makeReport {
path "wf-transcriptomes-report.html", emit: report
script:
String metadata_json = new JsonBuilder(metadata).toPrettyString().replaceAll("'", "'\\\\''")
def report_stats = stats instanceof java.util.Collection ? stats : (stats ? [stats] : [])
def report_pychopper = pychopper_dirs instanceof java.util.Collection ? pychopper_dirs : (pychopper_dirs ? [pychopper_dirs] : [])
def report_stats = (stats instanceof java.util.Collection ? stats : (stats ? [stats] : []))
.findAll { it.name != OPTIONAL_FILE.name }
String stats_args = report_stats ? "--stats ${report_stats.join(' ')}" : ""
String pychopper_args = report_pychopper.find { it.name != OPTIONAL_FILE.name } ? "--pychopper_dir pychopper" : ""
String de_args = de_files.name == OPTIONAL_FILE.name ? "" : "--de_dir de_analysis"
"""
echo '${metadata_json}' > metadata.json
workflow-glue report wf-transcriptomes-report.html \
--metadata metadata.json \
${stats_args} \
--alignment_stats_dir alignment_stats \
--cohort_dir cohort \
--samples_dir samples \
${pychopper_args} \
--sqanti_dir sqanti \
${de_args} \
--versions versions \
@ -179,7 +112,6 @@ def coerceBooleanParam(value) {
"version",
"igv",
"direct_rna",
"cdna_preprocess",
"de_analysis",
"analyse_unclassified",
"analyse_fail",
@ -211,17 +143,10 @@ workflow pipeline {
sample_sheet
ref_genome
ref_annotation
pychopper_dirs
main:
software_versions = getVersions()
workflow_params = getParams()
ingress_results = collectIngressResultsInDir(
reads.map { meta, sample_reads, stats ->
[meta, sample_reads, stats ?: OPTIONAL_FILE]
}
)
transcriptome = transcriptome_analysis(reads, ref_genome, ref_annotation, sample_sheet)
if (params.de_analysis) {
@ -236,50 +161,53 @@ workflow pipeline {
}
report_input = reads
.map { meta, sample_reads, stats ->
["all_samples", meta + [has_stats: stats as boolean], stats]
.collect(flat: false)
.map { rows ->
def metadata = rows.collect { meta, xam, xai, stats ->
meta + [has_stats: stats as boolean]
}
def stats_rows = rows.findAll { meta, xam, xai, stats -> stats != null }
[
metadata,
stats_rows ? stats_rows.collect { meta, xam, xai, stats -> stats } : [OPTIONAL_FILE]
]
}
.groupTuple()
.map { report_group, metas, stats ->
[metas, stats.findAll { it != null }]
}
alignment_stats = transcriptome.alignments
.map { meta, bam, bai, flagstat -> flagstat }
.collect()
sample_dirs_for_report = transcriptome.sample_dirs
.map { meta, sample_dir -> sample_dir }
.collect()
pychopper_dirs_for_report = pychopper_dirs
.map { meta, pychopper_dir -> pychopper_dir }
.ifEmpty(OPTIONAL_FILE)
.collect()
sqanti_dirs_for_report = transcriptome.joint_sqanti_dir
.concat(transcriptome.sample_sqanti_dirs.map { meta, sqanti_dir -> sqanti_dir })
.ifEmpty(OPTIONAL_FILE)
.collect()
// meta.src_xam is non-null if BAMs are "passed through"
generated_alignment_outputs = reads
.filter { meta, bam, bai, stats -> meta.src_xam == null }
.flatMap { meta, bam, bai, stats ->
def outdir = "cohort/alignments/${meta.alias}"
[
[bam, outdir],
[bai, outdir],
[stats.resolve("bamstats.flagstat.tsv"), outdir],
]
}
report = makeReport(
report_input,
software_versions,
workflow_params,
alignment_stats,
transcriptome.joint_dir,
sample_dirs_for_report,
pychopper_dirs_for_report,
sqanti_dirs_for_report,
de_dir,
workflow.manifest.version
)
results = Channel.empty()
.concat(ingress_results.out.map { [it, "ingress_results"] })
.concat(report.report.map { [it, null] })
.concat(workflow_params.map { [it, null] })
.concat(pychopper_dirs.map { meta, pychopper_dir -> [pychopper_dir, "ingress_results/${meta.alias}"] })
.concat(transcriptome.annotation_reference_summary.map { [it, "cohort/reference"] })
.concat(transcriptome.unstranded_annotation.map { [it, "cohort/reference"] })
.concat(transcriptome.joint_gtf.map { [it, "cohort"] })
@ -296,27 +224,28 @@ workflow pipeline {
.concat(transcriptome.sample_transcript_rds.map { meta, rds -> [rds, "samples/${meta.alias}"] })
.concat(transcriptome.sample_gene_rds.map { meta, rds -> [rds, "samples/${meta.alias}"] })
.concat(transcriptome.sample_metadata.map { meta, metadata -> [metadata, "samples/${meta.alias}"] })
.concat(transcriptome.alignments.map { meta, bam, bai, flagstat -> [bam, "cohort/alignments"] })
.concat(transcriptome.alignments.map { meta, bam, bai, flagstat -> [bai, "cohort/alignments"] })
.concat(transcriptome.alignments.map { meta, bam, bai, flagstat -> [flagstat, "cohort/alignments"] })
.concat(generated_alignment_outputs)
.concat(transcriptome.joint_sqanti_dir.map { [it, "cohort"] })
.concat(transcriptome.sample_sqanti_dirs.map { meta, sqanti_dir -> [sqanti_dir, "samples/${meta.alias}"] })
reference_basename = file(params.ref_genome).getName()
if (params.igv) {
results = results
.concat(transcriptome.reference.map { [it, "igv_reference"] })
.concat(transcriptome.reference_fai.map { [it, "igv_reference"] })
.concat(transcriptome.reference_gzi.map { [it, "igv_reference"] })
//todo ???
//results = results
// .concat(transcriptome.reference.map { [it, "igv_reference"] })
// .concat(transcriptome.reference_fai.map { [it, "igv_reference"] })
// .concat(transcriptome.reference_gzi.map { [it, "igv_reference"] })
results = Channel.empty()
igv_index_paths = transcriptome.reference_fai
.map { "igv_reference/${it.getName()}" }
.concat(transcriptome.reference_gzi.map { "igv_reference/${it.getName()}" })
//igv_index_paths = transcriptome.reference_fai
// .map { "igv_reference/${it.getName()}" }
// .concat(transcriptome.reference_gzi.map { "igv_reference/${it.getName()}" })
igv_index_paths = Channel.empty()
igv_alignment_paths = transcriptome.alignments
.map { meta, bam, bai, flagstat -> [
"cohort/alignments/${bam.getName()}",
"cohort/alignments/${bai.getName()}"
igv_alignment_paths = reads
.map { meta, bam, bai, stat -> [
meta.src_xam ?: "cohort/alignments/${meta.alias}/reads.bam",
meta.src_xai ?: "cohort/alignments/${meta.alias}/reads.bam.bai"
] }
.flatten()
@ -358,7 +287,7 @@ workflow {
throw new Exception("Provide exactly one of --fastq or --bam.")
}
if (!params.ref_genome) {
throw new Exception("Provide --ref_genome.")
throw new Exception("Provide --ref_genome.") //todo isnt this enforced in the schema?
}
if (!params.ref_annotation) {
throw new Exception("Provide --ref_annotation.")
@ -366,20 +295,20 @@ workflow {
if (!(params.transcriptome_mode in ["discover", "fixed_annotation"])) {
throw new Exception("--transcriptome_mode must be one of: discover, fixed_annotation.")
}
if (params.direct_rna && params.cdna_preprocess) {
throw new Exception("--cdna_preprocess cannot be used together with --direct_rna.")
}
if (params.de_analysis && !params.sample_sheet) {
throw new Exception("Provide --sample_sheet when running with --de_analysis.")
}
sample_sheet = params.sample_sheet ? file(params.sample_sheet, type: "file") : OPTIONAL_FILE
ref_genome = file(params.ref_genome, type: "file")
ref_annotation = file(params.ref_annotation, type: "file")
if (!ref_genome.exists()) {
throw new Exception("--ref_genome does not exist.")
}
prepared_reference = prepare_reference(
params.ref_genome, [
"output_cache": false,
"output_mmi": false,
])
ref_genome = prepared_reference.ref_tuple
if (!ref_annotation.exists()) {
throw new Exception("--ref_annotation does not exist.")
}
@ -387,51 +316,31 @@ workflow {
throw new Exception("--sample_sheet does not exist.")
}
def samples
def ingress_args = [
"minimap2_memory": ["31GB", "62GB"],
"minimap2_opts": params.direct_rna ? "-ax splice -uf -k14" : "-ax splice -uf",
"alignment_threads": 12,
"output_xam_fmt": "bam",
"sample": params.sample,
"sample_sheet": params.sample_sheet,
"analyse_unclassified": params.analyse_unclassified,
"analyse_fail": params.analyse_fail,
"fastcat_extra_args": "",
"required_sample_types": [],
]
if (params.fastq) {
samples = fastq_ingress([
"input": params.fastq,
"sample": params.sample,
"sample_sheet": params.sample_sheet,
"analyse_unclassified": params.analyse_unclassified,
"analyse_fail": params.analyse_fail,
"fastcat_extra_args": "",
"required_sample_types": [],
"fastq_chunk": params.fastq_chunk,
"per_read_stats": params.wf.per_read_stats,
"allow_multiple_basecall_models": params.wf.allow_multiple_basecall_models,
])
] + ingress_args, ref_genome)
} else {
samples = xam_ingress([
"input": params.bam,
"sample": params.sample,
"sample_sheet": params.sample_sheet,
"analyse_unclassified": params.analyse_unclassified,
"analyse_fail": params.analyse_fail,
"keep_unaligned": params.wf.keep_unaligned,
"return_fastq": params.wf.return_fastq,
"fastq_chunk": params.fastq_chunk,
"per_read_stats": params.wf.per_read_stats,
"allow_multiple_basecall_models": params.wf.allow_multiple_basecall_models,
])
] + ingress_args, ref_genome)
}
decorated_samples = samples
.map { meta, fname, stats -> [meta["group_key"], meta, fname, stats] }
.groupTuple()
.map { key, metas, fnames, statss ->
if (fnames[0] == null) {
fnames = null
}
[
metas[0] + ["group_index": metas.collect { it["group_index"] }],
fnames,
statss[0]
]
}
analysis_samples = decorated_samples
.filter { meta, sample_reads, stats ->
analysis_samples = samples
.filter { meta, xam, xai, stats ->
if (meta.n_seqs == 0) {
log.warn("Sample ${meta.alias} has no reads - excluded from transcriptome analysis.")
return false
@ -442,39 +351,10 @@ workflow {
throw new Exception("No samples with reads were available for transcriptome analysis.")
}
pychopper_results = Channel.empty()
processed_samples = analysis_samples
if (params.cdna_preprocess) {
grouped_samples = analysis_samples.branch { meta, sample_reads, stats ->
to_process: sample_reads != null
passthrough: sample_reads == null
}
preprocessed = preprocess_reads(
grouped_samples.to_process
.map { meta, sample_reads, stats -> [meta, sample_reads] }
)
processed_samples = grouped_samples.passthrough
.mix(
grouped_samples.to_process
.map { meta, sample_reads, stats -> [meta.alias, meta, stats] }
.join(preprocessed.full_len_reads)
.map { alias, meta, stats, full_length_reads ->
[meta, full_length_reads, stats]
}
)
pychopper_results = grouped_samples.to_process
.map { meta, sample_reads, stats -> [meta.alias, meta] }
.join(preprocessed.dir)
.map { alias, meta, pychopper_dir ->
[meta, pychopper_dir]
}
}
pipeline_run = pipeline(processed_samples, sample_sheet, ref_genome, ref_annotation, pychopper_results)
pipeline_run = pipeline(processed_samples, sample_sheet, ref_genome, ref_annotation)
publishResults(pipeline_run.results)
}

4
nextflow.config
View File

@ -13,8 +13,6 @@ params {
ref_annotation = null
transcriptome_mode = "discover"
direct_rna = false
cdna_preprocess = false
cdna_kit = "SQK-PCS109"
de_analysis = false
condition_column = "condition"
@ -32,8 +30,6 @@ params {
skip_sqanti = false
sqanti_skip_orf = true
sqanti_extra_args = ""
pychopper_backend = "edlib"
pychopper_opts = ""
aws_image_prefix = null
aws_queue = null

26
nextflow_schema.json
View File

@ -96,12 +96,6 @@
"type": "boolean",
"default": false,
"description": "Set this for direct RNA sequencing libraries."
},
"cdna_preprocess": {
"type": "boolean",
"default": false,
"description": "Apply the optional cDNA preprocessing stage before alignment.",
"help_text": "This is only relevant for cDNA libraries and must not be combined with `--direct_rna`."
}
},
"required": [
@ -199,26 +193,6 @@
"type": "number",
"description": "Optional bambu novel discovery rate override."
},
"cdna_kit": {
"type": "string",
"default": "SQK-PCS109",
"description": "ONT cDNA kit identifier used for pychopper preprocessing.",
"help_text": "Only used when `--cdna_preprocess` is enabled. The workflow derives the pychopper kit code from this value."
},
"pychopper_backend": {
"type": "string",
"default": "edlib",
"enum": [
"edlib",
"phmm"
],
"description": "Primer-detection backend to use for pychopper preprocessing."
},
"pychopper_opts": {
"type": "string",
"default": "",
"description": "Extra command-line options to pass to pychopper."
},
"skip_sqanti": {
"type": "boolean",
"default": false,

10
output_definition.json
View File

@ -41,7 +41,7 @@
"type": "per-sample"
},
"sample-bam": {
"filepath": "cohort/alignments/{{ alias }}.aligned.sorted.bam",
"filepath": "cohort/alignments/{{ alias }}/reads.bam",
"title": "Aligned BAM",
"description": "Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV.",
"mime-type": "application/gzip",
@ -49,7 +49,7 @@
"type": "per-sample"
},
"sample-bai": {
"filepath": "cohort/alignments/{{ alias }}.aligned.sorted.bam.bai",
"filepath": "cohort/alignments/{{ alias }}/reads.bam.bai",
"title": "Aligned BAM index",
"description": "Index for the aligned BAM.",
"mime-type": "application/octet-stream",
@ -57,10 +57,10 @@
"type": "per-sample"
},
"sample-flagstat": {
"filepath": "cohort/alignments/{{ alias }}.flagstat.txt",
"filepath": "cohort/alignments/{{ alias }}/bamstats.flagstat.tsv",
"title": "Alignment summary",
"description": "samtools flagstat output for the aligned BAM.",
"mime-type": "text/plain",
"description": "bamstats flagstat summary for the aligned BAM.",
"mime-type": "text/tab-separated-values",
"optional": false,
"type": "per-sample"
},

132
subworkflows/transcriptome.nf
View File

@ -9,7 +9,7 @@ process prepareAnnotationReference {
memory "6 GB"
input:
path ref_annotation
path ref_genome
tuple path(ref), path(ref_idx)
output:
stdout emit: warnings
path "annotation.gtf", emit: annotation
@ -20,62 +20,19 @@ process prepareAnnotationReference {
"""
workflow-glue prepare_annotation_reference \
--annotation "${ref_annotation}" \
--reference "${ref_genome}" \
--reference "${ref}" \
--out_dir prepared
mv prepared/* .
"""
}
process buildMinimapIndex {
label "wf_transcriptomes"
cpus 4
memory "32 GB"
input:
path reference
output:
tuple path("genome.mmi"), path(reference), emit: index
script:
String index_opts = params.minimap2_index_opts ?: ""
"""
minimap2 -t ${task.cpus} ${index_opts} -d genome.mmi "${reference}"
"""
}
process alignReads {
label "wf_transcriptomes"
cpus { params.threads ?: 4 }
memory "16 GB"
input:
tuple val(meta), path(reads), path(index), path(reference)
output:
tuple val(meta),
path("${meta.alias}.aligned.sorted.bam"),
path("${meta.alias}.aligned.sorted.bam.bai"),
path("${meta.alias}.flagstat.txt"),
emit: bam
script:
String preset = params.direct_rna ? "-ax splice -uf -k14" : "-ax splice -uf"
String extra = params.minimap2_opts ?: ""
def read_args = reads instanceof java.util.ArrayList ? reads.join(" ") : reads
int sort_threads = Math.max(task.cpus as int - 1, 1)
"""
minimap2 -t ${task.cpus} ${preset} ${extra} "${index}" ${read_args} \
| samtools sort --write-index -@ ${sort_threads} \
-o "${meta.alias}.aligned.sorted.bam##idx##${meta.alias}.aligned.sorted.bam.bai" -
samtools flagstat "${meta.alias}.aligned.sorted.bam" > "${meta.alias}.flagstat.txt"
"""
}
process runJointBambu {
label "wf_transcriptomes"
cpus { params.threads ?: 4 }
memory "32 GB"
input:
path bam_files, stageAs: "bams/*"
path bam_indexes, stageAs: "bams/*"
tuple val(aliases), path(bams, stageAs: "bams/??.bam"), path(bais, stageAs: "bams/??.bam.bai")
path sample_sheet
path annotation, stageAs: "annotation/*"
path reference, stageAs: "reference/*"
@ -88,11 +45,16 @@ process runJointBambu {
path "cohort/bambu_genes.rds", emit: gene_rds
path "cohort/transcript_metadata.tsv", emit: transcript_metadata
script:
def bam_list = bams instanceof Collection ? bams : [bams] // todo dont run joint on single sample anyway
def alias_list = aliases instanceof Collection ? aliases : [aliases]
String bams_arg = "--bams '${bam_list.join(",")}'"
String aliases_arg = "--aliases '${alias_list.join(",")}'"
String sample_sheet_arg = sample_sheet.name == OPTIONAL_FILE.name ? "" : "--sample_sheet ${sample_sheet}"
String ndr_arg = params.ndr != null ? "--ndr ${params.ndr}" : ""
"""
supeRglue bambu \
--bam_dir bams \
${bams_arg} \
${aliases_arg} \
${sample_sheet_arg} \
--annotation "${annotation}" \
--genome "${reference}" \
@ -109,7 +71,7 @@ process runPerSampleBambu {
cpus { params.threads ?: 4 }
memory "24 GB"
input:
tuple val(meta), path(bam), path(bai), path(flagstat)
tuple val(meta), path(bam), path(bai), path(stats)
path annotation, stageAs: "annotation/*"
path reference, stageAs: "reference/*"
output:
@ -121,11 +83,13 @@ process runPerSampleBambu {
tuple val(meta), path("${meta.alias}/bambu_genes.rds"), emit: gene_rds
tuple val(meta), path("${meta.alias}/transcript_metadata.tsv"), emit: transcript_metadata
script:
String bams_arg = "--bams '${bam.toString()}'"
String aliases_arg = "--aliases '${meta.alias}'"
String ndr_arg = params.ndr != null ? "--ndr ${params.ndr}" : ""
"""
supeRglue bambu \
--bam_path "${bam}" \
--sample_alias "${meta.alias}" \
${bams_arg} \
${aliases_arg} \
--annotation "${annotation}" \
--genome "${reference}" \
--transcriptome_mode "${params.transcriptome_mode}" \
@ -234,40 +198,9 @@ process runPerSampleSqanti {
}
process faidx {
label "wf_transcriptomes"
cpus 1
memory "4 GB"
input:
path ref
output:
path("${ref}.fai"), emit: fai
script:
"""
samtools faidx "${ref}"
"""
}
process gzFaidx {
label "wf_transcriptomes"
cpus 1
memory "4 GB"
errorStrategy "ignore"
input:
path ref
output:
tuple path("${ref}.fai"), path("${ref}.gzi"), emit: indexes
script:
"""
samtools faidx "${ref}"
"""
}
workflow transcriptome_analysis {
take:
reads
alignments
ref_genome
ref_annotation
sample_sheet
@ -281,22 +214,23 @@ workflow transcriptome_analysis {
analysis_annotation = prepared_reference_annotation.annotation
analysis_reference = prepared_reference_annotation.reference
genome_index = buildMinimapIndex(analysis_reference)
aligned = alignReads(
reads.map { meta, sample_reads, stats -> [meta, sample_reads] }
.combine(genome_index.index)
)
joint_bambu = runJointBambu(
aligned.bam.map { meta, bam, bai, flagstat -> bam }.collect(),
aligned.bam.map { meta, bam, bai, flagstat -> bai }.collect(),
alignments
| collect(flat: false)
| map { rows ->
// transform [meta, bam, bai] to [[alias1...aliasN], [bam1...bamN], [bai1...baiN]]
tuple(
rows.collect { it[0].alias },
rows.collect { it[1] },
rows.collect { it[2] }
)
},
sample_sheet,
analysis_annotation,
analysis_reference
)
sample_bambu = runPerSampleBambu(aligned.bam, analysis_annotation, analysis_reference)
sample_bambu = runPerSampleBambu(alignments, analysis_annotation, analysis_reference)
joint_fasta = buildCohortTranscriptomeFasta(joint_bambu.gtf, analysis_reference)
sample_fastas = buildSampleTranscriptomeFasta(sample_bambu.gtf, analysis_reference)
@ -311,22 +245,10 @@ workflow transcriptome_analysis {
sample_sqanti_dirs = sample_sqanti.dir
}
compressed_reference = params.ref_genome.toLowerCase().endsWith("gz")
if (compressed_reference) {
ref_indexes = gzFaidx(ref_genome)
reference_fai = ref_indexes.indexes.map { fai, gzi -> fai }
reference_gzi = ref_indexes.indexes.map { fai, gzi -> gzi }
} else {
ref_indexes = faidx(ref_genome)
reference_fai = ref_indexes.fai
reference_gzi = Channel.empty()
}
emit:
reference = ref_genome
annotation = ref_annotation
annotation_reference_summary = prepared_reference_annotation.summary
unstranded_annotation = prepared_reference_annotation.unstranded
alignments = aligned.bam
joint_dir = joint_bambu.dir
joint_gtf = joint_bambu.gtf
joint_fasta = joint_fasta.fasta
@ -345,6 +267,4 @@ workflow transcriptome_analysis {
sample_metadata = sample_bambu.transcript_metadata
joint_sqanti_dir = joint_sqanti_dir
sample_sqanti_dirs = sample_sqanti_dirs
reference_fai = reference_fai
reference_gzi = reference_gzi
}