Simplify reference handling through the workflow prepareAnnotation [CW-7232]

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"],
     )
 
 

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,9 +94,6 @@ 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 annotation
@@ -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(

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:

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(