From 3d7d49795edf26b2bb3c6bc3231e398cf6d4d010 Mon Sep 17 00:00:00 2001 From: Chris Wright Date: Sat, 16 May 2026 13:22:47 +0000 Subject: [PATCH] [CW-7235] cast ? to . in gffs --- .gitlab-ci.yml | 8 + .../prepare_annotation_reference.py | 168 +++++++++++------- bin/workflow_glue/tests/common/conftest.py | 9 + .../test_check_experiment_design.py | 0 .../tests/{ => common}/test_report.py | 0 .../tests/{ => common}/test_report_qc.py | 0 .../{ => common}/test_summarise_sqanti.py | 0 bin/workflow_glue/tests/wf/conftest.py | 9 + .../test_prepare_annotation_reference.py | 105 ++++++----- nextflow.config | 2 +- 10 files changed, 190 insertions(+), 111 deletions(-) create mode 100644 bin/workflow_glue/tests/common/conftest.py rename bin/workflow_glue/tests/{ => common}/test_check_experiment_design.py (100%) rename bin/workflow_glue/tests/{ => common}/test_report.py (100%) rename bin/workflow_glue/tests/{ => common}/test_report_qc.py (100%) rename bin/workflow_glue/tests/{ => common}/test_summarise_sqanti.py (100%) create mode 100644 bin/workflow_glue/tests/wf/conftest.py rename bin/workflow_glue/tests/{ => wf}/test_prepare_annotation_reference.py (82%) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 5b66e5a..929a88f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -13,10 +13,18 @@ variables: CI_FLAVOUR: "new" PYTEST_CONTAINER_NAME: "wf-common" PYTEST_CONTAINER_CONFIG_KEY: "common_sha" + PYTEST_TESTS_PATH: "bin/workflow_glue/tests/common" RTEST_CONTAINER_NAME: "wf-transcriptomes-core" RTEST_CONTAINER_CONFIG_KEY: "container_sha" WF_TEMPLATE_ENFORCEMENT_BRANCH: "CW-6552" +pytest_wfcontainer: + extends: pytest + variables: + PYTEST_CONTAINER_NAME: "wf-transcriptomes-core" + PYTEST_CONTAINER_CONFIG_KEY: "container_sha" + PYTEST_TESTS_PATH: "bin/workflow_glue/tests/wf" + macos-run: # Let's avoid those ARM64 runners for now tags: diff --git a/bin/workflow_glue/prepare_annotation_reference.py b/bin/workflow_glue/prepare_annotation_reference.py index b047344..fcbe288 100644 --- a/bin/workflow_glue/prepare_annotation_reference.py +++ b/bin/workflow_glue/prepare_annotation_reference.py @@ -165,17 +165,100 @@ class Annotation: MAX_UNSTRANDED_EXAMPLES = 20 - def __init__(self, input_path, work_dir, gffread="gffread"): + def __init__(self, input_path, work_dir): """Initialize and convert annotation to GTF.""" self.input_path = Path(input_path) self.work_dir = Path(work_dir) - self.is_gff = _is_gff(input_path) - self.gffread = gffread - + self.normalised_records = 0 # outputs - self.unfiltered_path = None - self.prepared_path = None - self._convert_to_gtf() + self.was_gff = _is_gff(self.input_path) + self.unfiltered_path = self.input_path # intermediate GTF before filtering + + if not self.was_gff: + if not self.input_path.exists() or self.input_path.stat().st_size < 1: + raise ValueError(f"Prepared annotation is empty: {self.input_path}") + else: + intermediate = self._decompress_gff_if_needed(self.input_path) + intermediate, self.normalised_records = ( + self._normalise_unknown_gff_strands(intermediate) + ) + + self.unfiltered_path = self.work_dir / "annotation_unfiltered.gtf" + self._run_gffread_conversion(intermediate, self.unfiltered_path) + + self.output_path = self.unfiltered_path # maybe mutated after filtering + + def _normalise_unknown_gff_strands(self, source_path, dest_path=None): + """Replace '?' strand values with '.' and return (path_used, count).""" + source_path = Path(source_path) + + if dest_path is None: + dest_path = self.work_dir / "annotation_input_sanitised.gff" + dest_path = Path(dest_path) + normalised = 0 + + with open(source_path, encoding="utf-8") as src, open( + dest_path, "w", encoding="utf-8" + ) as dst: + for line in src: + stripped = line.rstrip("\n") + if not stripped or stripped.startswith("#"): + dst.write(line) + continue + + fields = stripped.split("\t") + if len(fields) >= 7 and fields[6] == "?": + fields[6] = "." + normalised += 1 + dst.write("\t".join(fields) + "\n") + + if normalised == 0: + dest_path.unlink() + return source_path, 0 + + return dest_path, normalised + + def _decompress_gff_if_needed(self, source_path=None): + """Materialise gzipped GFF input to plain text; pass through if plain.""" + if source_path is None: + source_path = self.input_path + source_path = Path(source_path) + + if not _is_gzip(source_path): + return source_path + + # gffread does not support gzip input + name_without_gz = source_path.name[:-3] + suffix = Path(name_without_gz).suffix or ".gff" + decompressed = self.work_dir / f"annotation_input{suffix}" + with gzip.open(source_path, "rb") as src, open( + decompressed, "wb" + ) as dst: + shutil.copyfileobj(src, dst) + + return decompressed + + def _run_gffread_conversion(self, gff_path, gtf_path): + """Run gffread conversion from GFF/GFF3 to GTF.""" + result = subprocess.run( + ["gffread", "-T", str(gff_path), "-o", str(gtf_path)], + check=False, + capture_output=True, + text=True, + ) + if result.returncode != 0: + details = "\n".join( + part + for part in (result.stdout, result.stderr) + if part + ) + raise ValueError( + "Failed to convert annotation to GTF with gffread.\n" + f"Return code: {result.returncode}\n" + "Error message:" + details + ) + if not Path(gtf_path).exists() or Path(gtf_path).stat().st_size < 1: + raise ValueError(f"Prepared annotation is empty: {gtf_path}") def _open_gtf(self, path): """Open a GTF file, handling gzip transparently.""" @@ -219,58 +302,10 @@ class Annotation: ids.add(fields[0]) return sorted(ids) - def _convert_to_gtf(self): - """Convert GFF/GFF3 to GTF using gffread, or use input if already GTF.""" - if self.is_gff: - dest = self.work_dir / "annotation_unfiltered.gtf" - - if _is_gzip(self.input_path): - # decompress as gff doesn't support gzip input - name_without_gz = self.input_path.name[:-3] - suffix = Path(name_without_gz).suffix or ".gff" - intermediate = self.work_dir / f"annotation_input{suffix}" - with gzip.open(self.input_path, "rb") as src, open( - intermediate, "wb" - ) as dst: - shutil.copyfileobj(src, dst) - else: - intermediate = self.input_path - - # use gffread to convert to GTF - result = subprocess.run( - [self.gffread, "-T", str(intermediate), "-o", str(dest)], - check=False, - capture_output=True, - text=True, - ) - if result.returncode != 0: - details = "\n".join( - part - for part in (result.stdout, result.stderr) - if part - ) - raise ValueError( - "Failed to convert annotation to GTF with gffread.\n" - f"Return code: {result.returncode}\n" - "Error message:" + details - ) - - if not dest.exists() or dest.stat().st_size < 1: - raise ValueError(f"Prepared annotation is empty: {dest}") - - self.unfiltered_path = dest - else: - # already GTF - if not self.input_path.exists() or self.input_path.stat().st_size < 1: - raise ValueError(f"Prepared annotation is empty: {self.input_path}") - self.unfiltered_path = self.input_path - @property def seqnames(self): - """Get seqnames from prepared (stranded) annotation.""" - if self.prepared_path is None: - raise ValueError("Must call filter_stranded() before accessing seqnames") - return self._extract_seqnames(self.prepared_path) + """Get seqnames from the current annotation output path.""" + return self._extract_seqnames(self.output_path) def filter_stranded(self): """Filter to stranded records, return statistics.""" @@ -329,7 +364,7 @@ class Annotation: if stats["excluded_unstranded_records"] == 0 and unstranded_dest.exists(): unstranded_dest.unlink() - self.prepared_path = dest + self.output_path = dest stats["unstranded_path"] = ( str(unstranded_dest) if stats["excluded_unstranded_records"] > 0 else None ) @@ -434,7 +469,7 @@ def build_warnings(filter_stats, seqnames, ref_hints, ann_hints): return warnings -def prepare_annotation_reference(annotation, reference, out_dir, gffread="gffread"): +def prepare_annotation_reference(annotation, reference, out_dir): """Prepare annotation.gtf and reference.fasta in out_dir.""" out_dir = Path(out_dir) out_dir.mkdir(parents=True, exist_ok=False) @@ -443,7 +478,7 @@ def prepare_annotation_reference(annotation, reference, out_dir, gffread="gffrea ref = PreparedReference(reference, out_dir) # prepare annotation (convert to GTF if needed, then filter to stranded) - ann = Annotation(annotation, out_dir, gffread=gffread) + ann = Annotation(annotation, out_dir) filter_stats = ann.filter_stranded() # check seqname overlap between prepared annotation and reference @@ -464,14 +499,21 @@ def prepare_annotation_reference(annotation, reference, out_dir, gffread="gffrea ref_hints = ref.detect_hints() ann_hints = ann.detect_hints() warnings = build_warnings(filter_stats, seqnames, ref_hints, ann_hints) + if ann.normalised_records: + warnings.append( + "Warning: Replaced " + f"{ann.normalised_records} " + "GFF records with unknown strand '?' to '.' before gffread conversion." + ) # build a summary summary = { "annotation": { "input": str(ann.input_path), - "prepared": str(ann.prepared_path), + "prepared": str(ann.output_path), "unfiltered": str(ann.unfiltered_path), - "was_gff": ann.is_gff, + "was_gff": ann.was_gff, + "normalised_unknown_strand_records": ann.normalised_records, **filter_stats, }, "reference": { @@ -505,7 +547,6 @@ def main(args): args.annotation, args.reference, args.out_dir, - gffread=args.gffread, ) for warning in summary["warnings"]: @@ -530,5 +571,4 @@ def argparser(): required=True, help="Output directory for prepared annotation/reference files.", ) - parser.add_argument("--gffread", default="gffread", help="gffread executable.") return parser diff --git a/bin/workflow_glue/tests/common/conftest.py b/bin/workflow_glue/tests/common/conftest.py new file mode 100644 index 0000000..6b04a68 --- /dev/null +++ b/bin/workflow_glue/tests/common/conftest.py @@ -0,0 +1,9 @@ +#!/usr/bin/env python +"""Import path setup for common workflow_glue tests.""" + +from pathlib import Path +import sys + + +# Add /host/bin so `import workflow_glue` resolves in CI containers. +sys.path.insert(0, str(Path(__file__).resolve().parents[3])) diff --git a/bin/workflow_glue/tests/test_check_experiment_design.py b/bin/workflow_glue/tests/common/test_check_experiment_design.py similarity index 100% rename from bin/workflow_glue/tests/test_check_experiment_design.py rename to bin/workflow_glue/tests/common/test_check_experiment_design.py diff --git a/bin/workflow_glue/tests/test_report.py b/bin/workflow_glue/tests/common/test_report.py similarity index 100% rename from bin/workflow_glue/tests/test_report.py rename to bin/workflow_glue/tests/common/test_report.py diff --git a/bin/workflow_glue/tests/test_report_qc.py b/bin/workflow_glue/tests/common/test_report_qc.py similarity index 100% rename from bin/workflow_glue/tests/test_report_qc.py rename to bin/workflow_glue/tests/common/test_report_qc.py diff --git a/bin/workflow_glue/tests/test_summarise_sqanti.py b/bin/workflow_glue/tests/common/test_summarise_sqanti.py similarity index 100% rename from bin/workflow_glue/tests/test_summarise_sqanti.py rename to bin/workflow_glue/tests/common/test_summarise_sqanti.py diff --git a/bin/workflow_glue/tests/wf/conftest.py b/bin/workflow_glue/tests/wf/conftest.py new file mode 100644 index 0000000..f45ec8a --- /dev/null +++ b/bin/workflow_glue/tests/wf/conftest.py @@ -0,0 +1,9 @@ +#!/usr/bin/env python +"""Import path setup for wf-container workflow_glue tests.""" + +from pathlib import Path +import sys + + +# Add /host/bin so `import workflow_glue` resolves in CI containers. +sys.path.insert(0, str(Path(__file__).resolve().parents[3])) diff --git a/bin/workflow_glue/tests/test_prepare_annotation_reference.py b/bin/workflow_glue/tests/wf/test_prepare_annotation_reference.py similarity index 82% rename from bin/workflow_glue/tests/test_prepare_annotation_reference.py rename to bin/workflow_glue/tests/wf/test_prepare_annotation_reference.py index 7b29f69..52f559f 100644 --- a/bin/workflow_glue/tests/test_prepare_annotation_reference.py +++ b/bin/workflow_glue/tests/wf/test_prepare_annotation_reference.py @@ -2,7 +2,6 @@ import gzip from pathlib import Path -from types import SimpleNamespace import pytest from workflow_glue import get_components @@ -31,14 +30,12 @@ chr1\tsim\texon\t1\t5\t.\t+\t.\tParent=transcript1 chr1\tsim\texon\t6\t10\t.\t+\t.\tParent=transcript1 """ -GFF3_AS_GTF = ( - 'chr1\tsim\ttranscript\t1\t10\t.\t+\t.\tgene_id "gene1"; ' - 'transcript_id "transcript1";\n' - 'chr1\tsim\texon\t1\t5\t.\t+\t.\tgene_id "gene1"; ' - 'transcript_id "transcript1";\n' - 'chr1\tsim\texon\t6\t10\t.\t+\t.\tgene_id "gene1"; ' - 'transcript_id "transcript1";\n' -) +GFF3_INPUT_WITH_UNKNOWN_STRAND = """##gff-version 3 +chr1\tsim\tgene\t1\t10\t.\t?\t.\tID=gene1 +chr1\tsim\tmRNA\t1\t10\t.\t+\t.\tID=transcript1;Parent=gene1 +chr1\tsim\texon\t1\t5\t.\t+\t.\tParent=transcript1 +chr1\tsim\texon\t6\t10\t.\t+\t.\tParent=transcript1 +""" NCBI_REFERENCE = """>chr1 GRCh38 RefSeq primary assembly AAAA @@ -72,27 +69,27 @@ MOUSE_GTF = ( @pytest.mark.parametrize( - "ref_data,ref_name,annot_data,annot_name,gff_output", + "ref_data,ref_name,annot_data,annot_name", [ (SIMPLE_REFERENCE, "reference.fa", SIMPLE_GTF, - "annotation.gtf", None), + "annotation.gtf"), (SIMPLE_REFERENCE, "reference.fa.gz", SIMPLE_GTF, - "annotation.gtf.gz", None), + "annotation.gtf.gz"), (SIMPLE_REFERENCE, "reference.fa", GFF3_INPUT, - "annotation.gff3.gz", GFF3_AS_GTF), + "annotation.gff3.gz"), (NCBI_REFERENCE, "reference.fna.gz", NCBI_GTF, - "annotation.gtf.gz", None), + "annotation.gtf.gz"), (SIMPLE_REFERENCE + ">chrExtra\nTTTT\n", "reference.fa", - MIXED_STRANDED_GTF, "annotation.gtf", None), + MIXED_STRANDED_GTF, "annotation.gtf"), (MOUSE_REFERENCE, "GRCm39.genome.fa", MOUSE_GTF, - "gencode.vM33.annotation.gtf", None), + "gencode.vM33.annotation.gtf"), ], ids=[ "simple_gtf", "gzipped_gtf", "gff3_input", "ncbi_format", "mixed_stranded", "mouse_gencode"], ) def test_prepare_all_formats_produce_valid_outputs( - ref_data, ref_name, annot_data, annot_name, gff_output, tmp_path, monkeypatch + ref_data, ref_name, annot_data, annot_name, tmp_path ): """All supported input formats produce annotation.gtf and reference.fasta.""" # write reference @@ -109,14 +106,6 @@ def test_prepare_all_formats_produce_valid_outputs( else: _write(annot_path, annot_data) - # mock gffread (only invoked for GFF3 inputs) - if gff_output: - def fake_run(command, check, capture_output, text): - output_idx = command.index("-o") + 1 - Path(command[output_idx]).write_text(gff_output, encoding="utf-8") - return SimpleNamespace(returncode=0, stdout="", stderr="") - monkeypatch.setattr(prepare_annotation_reference.subprocess, "run", fake_run) - out_dir = tmp_path / "prepared" summary = prepare_annotation_reference.prepare_annotation_reference( annot_path, ref_path, out_dir @@ -172,34 +161,59 @@ def test_ncbi_format_sanitises_gene_ids(tmp_path): ) -def test_gff3_conversion_via_gffread(tmp_path, monkeypatch): +def test_gff3_conversion_via_gffread(tmp_path): """GFF3 inputs invoke gffread and convert to GTF format.""" reference = _write(tmp_path / "reference.fa", SIMPLE_REFERENCE) annotation = _write(tmp_path / "annotation.gff3", GFF3_INPUT) - captured = {} - - def fake_run(command, check, capture_output, text): - captured["command"] = command - output_idx = command.index("-o") + 1 - Path(command[output_idx]).write_text(GFF3_AS_GTF, encoding="utf-8") - return SimpleNamespace(returncode=0, stdout="", stderr="") - monkeypatch.setattr(prepare_annotation_reference.subprocess, "run", fake_run) summary = prepare_annotation_reference.prepare_annotation_reference( annotation, reference, tmp_path / "prepared", - gffread="custom-gffread", ) assert summary["annotation"]["was_gff"] is True - assert captured["command"][0] == "custom-gffread" - assert "-T" in captured["command"] prepared_text = Path(summary["annotation"]["prepared"]).read_text(encoding="utf-8") assert 'gene_id "gene1"' in prepared_text assert 'transcript_id "transcript1"' in prepared_text - assert prepared_text.count("exon") == 2 + assert "\texon\t" in prepared_text + + +def test_normalise_unknown_gff_strands_rewrites_question_mark(tmp_path): + """Unknown GFF strand '?' should be rewritten to '.'.""" + input_path = _write(tmp_path / "annotation.gff3", GFF3_INPUT_WITH_UNKNOWN_STRAND) + output_path = tmp_path / "annotation_sanitised.gff" + init_annotation = _write(tmp_path / "annotation_init.gtf", SIMPLE_GTF) + ann = prepare_annotation_reference.Annotation( + init_annotation, + tmp_path, + ) + + used_path, normalised = ann._normalise_unknown_gff_strands(input_path, output_path) + + assert used_path == output_path + assert normalised == 1 + output_text = output_path.read_text(encoding="utf-8") + assert "\t?\t" not in output_text + assert "\t.\t" in output_text + + +def test_gff3_unknown_strand_with_real_gffread(tmp_path): + """Integration: real gffread accepts input after '?' strand normalisation.""" + reference = _write(tmp_path / "reference.fa", SIMPLE_REFERENCE) + annotation = _write(tmp_path / "annotation.gff3", GFF3_INPUT_WITH_UNKNOWN_STRAND) + + summary = prepare_annotation_reference.prepare_annotation_reference( + annotation, + reference, + tmp_path / "prepared", + ) + + prepared_text = Path(summary["annotation"]["prepared"]).read_text(encoding="utf-8") + assert summary["annotation"]["normalised_unknown_strand_records"] == 1 + assert "\t?\t" not in prepared_text + assert summary["annotation"]["kept_records"] > 0 def test_unstranded_records_filtered_and_saved_separately(tmp_path): @@ -250,16 +264,15 @@ def test_seqname_warnings_reflect_filtered_annotation(tmp_path): ) -def test_gffread_failure_raises_error(tmp_path, monkeypatch): - """Failed GFF3 conversion should abort before downstream processing.""" +def test_invalid_gff_for_conversion_raises_error(tmp_path): + """Malformed GFF input should fail during gffread conversion.""" reference = _write(tmp_path / "reference.fa", SIMPLE_REFERENCE) - annotation = _write(tmp_path / "annotation.gff3", GFF3_INPUT) + annotation = _write(tmp_path / "annotation.gff3", "chr1\tsim\tgene\t1\t10\n") - def fake_run(command, check, capture_output, text): - return SimpleNamespace(returncode=1, stdout="", stderr="gffread error") - monkeypatch.setattr(prepare_annotation_reference.subprocess, "run", fake_run) - - with pytest.raises(ValueError, match="Failed to convert annotation"): + with pytest.raises( + ValueError, + match="(Failed to convert annotation to GTF|Prepared annotation is empty)", + ): prepare_annotation_reference.prepare_annotation_reference( annotation, reference, diff --git a/nextflow.config b/nextflow.config index 276c34a..88c5500 100644 --- a/nextflow.config +++ b/nextflow.config @@ -56,7 +56,7 @@ params { "--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv'", ] common_sha = "sha21d552f9910c575766e5d465fcb7b52fefda4b79" - container_sha = "shaff0012055c9e1e71caf5b7857d717d17a3465a77" + container_sha = "sha02e44f706d88fa29d8344b78479f187db7eec4ec" pychopper_sha = "shaaaf20a5a0e76f9e18bad21af639a6b69e4a31a2f" sqanti_sha = "sha5bd775836492699e2537ebf846098eb117191d87" agent = null