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

Merge branch 'cw-7235' into 'dev'

Browse Source 
[CW-7235] cast ? to . in gffs

See merge request epi2melabs/workflows/wf-transcriptomes!261
This commit is contained in:
Chris Wright 2026-05-16 13:22:47 +00:00
commit e2166fe062
10 changed files with 190 additions and 111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
.gitlab-ci.yml
View File

@ -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:

168
bin/workflow_glue/prepare_annotation_reference.py
View File

@ -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

9
bin/workflow_glue/tests/common/conftest.py Normal file
View File

@ -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]))

0
bin/workflow_glue/tests/test_check_experiment_design.py → bin/workflow_glue/tests/common/test_check_experiment_design.py
View File

0
bin/workflow_glue/tests/test_report.py → bin/workflow_glue/tests/common/test_report.py
View File

0
bin/workflow_glue/tests/test_report_qc.py → bin/workflow_glue/tests/common/test_report_qc.py
View File

0
bin/workflow_glue/tests/test_summarise_sqanti.py → bin/workflow_glue/tests/common/test_summarise_sqanti.py
View File

9
bin/workflow_glue/tests/wf/conftest.py Normal file
View File

@ -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]))

105
bin/workflow_glue/tests/test_prepare_annotation_reference.py → bin/workflow_glue/tests/wf/test_prepare_annotation_reference.py
View File

@ -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,

2
nextflow.config
View File

@ -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