From dc0ebf3a52209d21bf9da24995c65ec40fb6d655 Mon Sep 17 00:00:00 2001 From: Natalia Garcia Date: Wed, 20 May 2026 12:49:08 +0000 Subject: [PATCH] Simplify reference handling through the workflow prepareAnnotation [CW-7232] --- .../prepare_annotation_reference.py | 27 ++++--------------- .../wf/test_prepare_annotation_reference.py | 17 +++++------- modules/local/bambu_chunked.nf | 4 +-- subworkflows/transcriptome.nf | 12 ++++----- 4 files changed, 19 insertions(+), 41 deletions(-) diff --git a/bin/workflow_glue/prepare_annotation_reference.py b/bin/workflow_glue/prepare_annotation_reference.py index fcbe288..f0259a7 100644 --- a/bin/workflow_glue/prepare_annotation_reference.py +++ b/bin/workflow_glue/prepare_annotation_reference.py @@ -51,33 +51,17 @@ def _find_patterns(text, patterns): class PreparedReference: """A reference genome prepared for analysis.""" - def __init__(self, input_path, work_dir): + def __init__(self, input_path): """Initialize and prepare reference genome.""" self.input_path = Path(input_path) - self.work_dir = Path(work_dir) - self.output_path = self.work_dir / "reference.fasta" self._seqnames = None - self._prepare() - - def _prepare(self): - """Prepare reference (decompress if gzipped).""" - if _is_gzip(self.input_path): - with gzip.open(self.input_path, "rb") as src, open( - self.output_path, "wb" - ) as dst: - shutil.copyfileobj(src, dst) - else: - shutil.copyfile(self.input_path, self.output_path) - - if not self.output_path.exists() or self.output_path.stat().st_size < 1: - raise ValueError(f"Prepared file is empty: {self.output_path}") @property def seqnames(self): """Extract and cache sorted FASTA sequence names.""" if self._seqnames is None: ids = set() - with open(self.output_path, encoding="utf-8") as f: + with open(self.input_path, encoding="utf-8") as f: for line in f: if line.startswith(">"): ids.add(line[1:].split()[0]) @@ -94,7 +78,7 @@ class PreparedReference: providers.update(_find_patterns(self.input_path.name, PROVIDER_PATTERNS)) # check headers - with open(self.output_path, encoding="utf-8") as f: + with open(self.input_path, encoding="utf-8") as f: for line in f: if line.startswith(">"): builds.update(_find_patterns(line, BUILD_PATTERNS)) @@ -475,7 +459,7 @@ def prepare_annotation_reference(annotation, reference, out_dir): out_dir.mkdir(parents=True, exist_ok=False) # prepare reference genome - ref = PreparedReference(reference, out_dir) + ref = PreparedReference(reference) # prepare annotation (convert to GTF if needed, then filter to stranded) ann = Annotation(annotation, out_dir) @@ -518,7 +502,6 @@ def prepare_annotation_reference(annotation, reference, out_dir): }, "reference": { "input": str(ref.input_path), - "prepared": str(ref.output_path), }, "seqnames": seqnames, "analysis_seqnames": seqnames, @@ -555,7 +538,7 @@ def main(args): logger.info( "Prepared annotation %s and reference %s.", summary["annotation"]["prepared"], - summary["reference"]["prepared"], + summary["reference"]["input"], ) diff --git a/bin/workflow_glue/tests/wf/test_prepare_annotation_reference.py b/bin/workflow_glue/tests/wf/test_prepare_annotation_reference.py index 52f559f..8c12367 100644 --- a/bin/workflow_glue/tests/wf/test_prepare_annotation_reference.py +++ b/bin/workflow_glue/tests/wf/test_prepare_annotation_reference.py @@ -73,11 +73,11 @@ MOUSE_GTF = ( [ (SIMPLE_REFERENCE, "reference.fa", SIMPLE_GTF, "annotation.gtf"), - (SIMPLE_REFERENCE, "reference.fa.gz", SIMPLE_GTF, + (SIMPLE_REFERENCE, "reference.fa", SIMPLE_GTF, "annotation.gtf.gz"), (SIMPLE_REFERENCE, "reference.fa", GFF3_INPUT, "annotation.gff3.gz"), - (NCBI_REFERENCE, "reference.fna.gz", NCBI_GTF, + (NCBI_REFERENCE, "reference.fna", NCBI_GTF, "annotation.gtf.gz"), (SIMPLE_REFERENCE + ">chrExtra\nTTTT\n", "reference.fa", MIXED_STRANDED_GTF, "annotation.gtf"), @@ -94,10 +94,7 @@ def test_prepare_all_formats_produce_valid_outputs( """All supported input formats produce annotation.gtf and reference.fasta.""" # write reference ref_path = tmp_path / ref_name - if ref_name.endswith(".gz"): - _write_gzip(ref_path, ref_data) - else: - _write(ref_path, ref_data) + _write(ref_path, ref_data) # write annotation annot_path = tmp_path / annot_name @@ -113,15 +110,15 @@ def test_prepare_all_formats_produce_valid_outputs( # output files must exist assert (out_dir / "annotation.gtf").exists() - assert (out_dir / "reference.fasta").exists() + assert (ref_path).exists() assert (out_dir / "annotation_reference_summary.json").exists() # paths returned in summary must match assert summary["annotation"]["prepared"] == str(out_dir / "annotation.gtf") - assert summary["reference"]["prepared"] == str(out_dir / "reference.fasta") + assert summary["reference"]["input"] == str(ref_path) # outputs must be valid for downstream tools (Bambu, SQANTI) - reference = Path(summary["reference"]["prepared"]) + reference = Path(summary["reference"]["input"]) annotation = Path(summary["annotation"]["prepared"]) assert reference.read_text(encoding="utf-8").startswith(">") assert summary["analysis_seqnames"]["has_overlap"] is True @@ -142,7 +139,7 @@ def test_prepare_all_formats_produce_valid_outputs( def test_ncbi_format_sanitises_gene_ids(tmp_path): """NCBI gene_id='transcript_id' should be sanitised to actual transcript_id.""" - reference = _write_gzip(tmp_path / "reference.fna.gz", NCBI_REFERENCE) + reference = _write(tmp_path / "reference.fna", NCBI_REFERENCE) annotation = _write_gzip(tmp_path / "annotation.gtf.gz", NCBI_GTF) summary = prepare_annotation_reference.prepare_annotation_reference( diff --git a/modules/local/bambu_chunked.nf b/modules/local/bambu_chunked.nf index be8a586..7d009db 100644 --- a/modules/local/bambu_chunked.nf +++ b/modules/local/bambu_chunked.nf @@ -14,7 +14,7 @@ process bambuDiscover { input: tuple val(meta), val(aliases), path(bams, stageAs: "bams/??.bam"), path(bais, stageAs: "bams/??.bam.bai"), path(sample_sheet) path annotation, stageAs: "annotation/*" - path reference, stageAs: "reference/*" + tuple path(reference, stageAs: "reference/reference.fa"), path(ref_fai, stageAs: "reference/reference.fai") output: tuple val(meta), path("discover"), emit: dir script: @@ -51,7 +51,7 @@ process bambuQuant { errorStrategy 'retry' input: tuple val(meta), val(chunk_id), val(annotation_tx_count), path(chunk_rds), path(discovered_annotation) - path reference, stageAs: "reference/*" + tuple path(reference, stageAs: "reference/reference.fa"), path(ref_fai, stageAs: "reference/reference.fai") output: tuple val(meta), val(chunk_id), path("${chunk_id}"), emit: dir script: diff --git a/subworkflows/transcriptome.nf b/subworkflows/transcriptome.nf index a8d4be7..6e99251 100644 --- a/subworkflows/transcriptome.nf +++ b/subworkflows/transcriptome.nf @@ -91,7 +91,6 @@ process prepareAnnotationReference { output: stdout emit: warnings path "annotation.gtf", emit: annotation - path "reference.fasta", emit: reference path "annotation_reference_summary.json", emit: summary path "unstranded_annotation.gtf", optional: true, emit: unstranded script: @@ -111,7 +110,7 @@ process buildCohortTranscriptomeFasta { memory "4 GB" input: path "transcripts.gtf" - path reference + tuple path(reference), path(ref_idx) output: path "cohort.transcriptome.fa", emit: fasta script: @@ -127,7 +126,7 @@ process buildSampleTranscriptomeFasta { memory "4 GB" input: tuple val(meta), path("transcripts.gtf") - path reference + tuple path(reference), path(ref_idx) output: tuple val(meta), path("${meta.alias}.transcriptome.fa"), emit: fasta script: @@ -144,7 +143,7 @@ process runJointSqanti { input: path gtf path annotation, stageAs: "annotation/*" - path reference, stageAs: "reference/*" + tuple path(reference, stageAs: "reference/reference.fa"), path(ref_fai, stageAs: "reference/reference.fai") output: path "sqanti_cohort", emit: dir path "sqanti_cohort/classification_summary.tsv", emit: summary @@ -177,7 +176,7 @@ process runPerSampleSqanti { input: tuple val(meta), path(gtf) path annotation, stageAs: "annotation/*" - path reference, stageAs: "reference/*" + tuple path(reference, stageAs: "reference/reference.fa"), path(ref_fai, stageAs: "reference/reference.fai") output: tuple val(meta), path("${meta.alias}_sqanti"), emit: dir tuple val(meta), path("${meta.alias}_sqanti/classification_summary.tsv"), emit: summary @@ -217,8 +216,7 @@ workflow transcriptome_analysis { } } analysis_annotation = prepared_reference_annotation.annotation.first() - analysis_reference = prepared_reference_annotation.reference.first() - + analysis_reference = ref_genome.first() joint_meta = [alias: "cohort"] joint_discover = runJointBambuDiscover(