diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d3e38c9..f920a93 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -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" diff --git a/README.md b/README.md index 507b086..5ecdc34 100644 --- a/README.md +++ b/README.md @@ -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//_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//`. -### 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//` for sample-specific models. If `--cdna_preprocess` is enabled, -look under `ingress_results//` for the `pychopper` outputs associated -with each sample. +`samples//` 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 diff --git a/bin/workflow_glue/models/common.py b/bin/workflow_glue/models/common.py index bef72fc..1028e8d 100644 --- a/bin/workflow_glue/models/common.py +++ b/bin/workflow_glue/models/common.py @@ -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) diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index 0651ab5..f13f1cb 100644 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -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, diff --git a/bin/workflow_glue/tests/test_report.py b/bin/workflow_glue/tests/test_report.py index 92f2e18..5dfa798 100644 --- a/bin/workflow_glue/tests/test_report.py +++ b/bin/workflow_glue/tests/test_report.py @@ -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, diff --git a/bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py b/bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py index ff7f3b4..b1cbabc 100644 --- a/bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py +++ b/bin/workflow_glue/wfg_helpers/check_bam_headers_in_dir.py @@ -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 diff --git a/bin/workflow_glue_r/R/bambu.R b/bin/workflow_glue_r/R/bambu.R index 251df83..3555dc9 100644 --- a/bin/workflow_glue_r/R/bambu.R +++ b/bin/workflow_glue_r/R/bambu.R @@ -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) diff --git a/bin/workflow_glue_r/R/de_analysis.R b/bin/workflow_glue_r/R/de_analysis.R index a3b9c97..e5ce5d7 100644 --- a/bin/workflow_glue_r/R/de_analysis.R +++ b/bin/workflow_glue_r/R/de_analysis.R @@ -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 ) }), diff --git a/bin/workflow_glue_r/tests/testthat/test_bambu.R b/bin/workflow_glue_r/tests/testthat/test_bambu.R index d803a24..2bd1140 100644 --- a/bin/workflow_glue_r/tests/testthat/test_bambu.R +++ b/bin/workflow_glue_r/tests/testthat/test_bambu.R @@ -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, diff --git a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R b/bin/workflow_glue_r/tests/testthat/test_de_analysis.R index de13f5b..d58dddd 100644 --- a/bin/workflow_glue_r/tests/testthat/test_de_analysis.R +++ b/bin/workflow_glue_r/tests/testthat/test_de_analysis.R @@ -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)) }) diff --git a/docs/02_introduction.md b/docs/02_introduction.md index f7fc250..72c427c 100644 --- a/docs/02_introduction.md +++ b/docs/02_introduction.md @@ -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` diff --git a/docs/06_pipeline_overview.md b/docs/06_pipeline_overview.md index 9cfc7a6..669ddaa 100644 --- a/docs/06_pipeline_overview.md +++ b/docs/06_pipeline_overview.md @@ -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//_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//`. -### 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: diff --git a/docs/07_input_parameters.md b/docs/07_input_parameters.md index 5425458..e1c877c 100644 --- a/docs/07_input_parameters.md +++ b/docs/07_input_parameters.md @@ -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. | | | diff --git a/docs/08_outputs.md b/docs/08_outputs.md index 3729e61..e7ec9ae 100644 --- a/docs/08_outputs.md +++ b/docs/08_outputs.md @@ -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 | diff --git a/docs/10_troubleshooting.md b/docs/10_troubleshooting.md index 17af1e5..6d0d996 100644 --- a/docs/10_troubleshooting.md +++ b/docs/10_troubleshooting.md @@ -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//` for sample-specific models. If `--cdna_preprocess` is enabled, -look under `ingress_results//` for the `pychopper` outputs associated -with each sample. +`samples//` 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 diff --git a/lib/common.nf b/lib/common.nf index 556d249..db6055c 100644 --- a/lib/common.nf +++ b/lib/common.nf @@ -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 ; + + """ +} diff --git a/lib/ingress.nf b/lib/ingress.nf index b6d3157..48150f4 100644 --- a/lib/ingress.nf +++ b/lib/ingress.nf @@ -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.. - 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.. + 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, diff --git a/lib/reference.nf b/lib/reference.nf new file mode 100644 index 0000000..b5c9235 --- /dev/null +++ b/lib/reference.nf @@ -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 +} diff --git a/main.nf b/main.nf index 676b9f7..504a214 100644 --- a/main.nf +++ b/main.nf @@ -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) } diff --git a/nextflow.config b/nextflow.config index 7fb2bcd..276c34a 100644 --- a/nextflow.config +++ b/nextflow.config @@ -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 diff --git a/nextflow_schema.json b/nextflow_schema.json index 3a67ed2..84f09c1 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -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, diff --git a/output_definition.json b/output_definition.json index f57a98e..31b3120 100644 --- a/output_definition.json +++ b/output_definition.json @@ -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" }, diff --git a/subworkflows/transcriptome.nf b/subworkflows/transcriptome.nf index d8d8fed..7d7f738 100644 --- a/subworkflows/transcriptome.nf +++ b/subworkflows/transcriptome.nf @@ -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 }