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

Improve messaging for seqnames with no transcript records [CW-7391]

Browse Source
This commit is contained in:
Sam Nicholls 2026-07-13 16:00:41 +00:00
parent c09fa66ffc
commit 4613584c69
6 changed files with 218 additions and 60 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CHANGELOG.md
View File

@ -7,6 +7,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
## [Unreleased]
### Added
- Annotation preparation summary JSON and output report include a count of non-transcript records and associated seqnames that are pruned when converting a GFF to GTF. Seqnames that do not appear in the converted GTF but were present in the input GFF are now classified as "only in GFF" rather than "only in reference".
### Fixed
- Incorrect strand assignment when mapping cDNA reads is fixed by using minimap2 `-ub` instead of `-uf`.
- Remove mention of analysis_group and type from `README.md` as they are not relevant for this workflow.

83
bin/workflow_glue/prepare_annotation_reference.py
View File

@ -154,6 +154,10 @@ class Annotation:
self.input_path = Path(input_path)
self.work_dir = Path(work_dir)
self.normalised_records = 0
self.gff_convert_input_count = None
self.gff_convert_output_count = None
self.gff_convert_pruned_count = None
self.only_in_gff = []
# outputs
self.was_gff = _is_gff(self.input_path)
self.unfiltered_path = self.input_path # intermediate GTF before filtering
@ -166,9 +170,23 @@ class Annotation:
intermediate, self.normalised_records = (
self._normalise_unknown_gff_strands(intermediate)
)
self.gff_convert_input_count = self._count_records(intermediate)
gff_input_seqnames = self._extract_seqnames(intermediate)
self.unfiltered_path = self.work_dir / "annotation_unfiltered.gtf"
self._run_gffread_conversion(intermediate, self.unfiltered_path)
self.gff_convert_output_count = self._count_records(
self.unfiltered_path
)
gff_converted_seqnames = self._extract_seqnames(
self.unfiltered_path
)
self.gff_convert_pruned_count = (
self.gff_convert_input_count - self.gff_convert_output_count
)
self.only_in_gff = sorted(
set(gff_input_seqnames) - set(gff_converted_seqnames)
)
self.output_path = self.unfiltered_path # maybe mutated after filtering
@ -286,6 +304,16 @@ class Annotation:
ids.add(fields[0])
return sorted(ids)
def _count_records(self, path):
"""Count non-empty, non-comment annotation records."""
records = 0
with self._open_gtf(path) as f:
for line in f:
stripped = line.strip()
if stripped and not stripped.startswith("#"):
records += 1
return records
@property
def seqnames(self):
"""Get seqnames from the current annotation output path."""
@ -374,13 +402,18 @@ class Annotation:
}
def validate_seqname_overlap(annotation_seqnames, reference_seqnames):
def validate_seqname_overlap(
annotation_seqnames,
reference_seqnames,
gff_seqnames=None,
):
"""Return overlap statistics for annotation/reference seqnames."""
ann_set = set(annotation_seqnames)
ref_set = set(reference_seqnames)
gff_set = set(gff_seqnames or [])
matches = sorted(ann_set & ref_set)
only_in_annotation = sorted(ann_set - ref_set)
only_in_reference = sorted(ref_set - ann_set)
only_in_reference = sorted(ref_set - ann_set - gff_set)
return {
"has_overlap": bool(matches),
@ -390,10 +423,20 @@ def validate_seqname_overlap(annotation_seqnames, reference_seqnames):
}
def build_warnings(filter_stats, seqnames, ref_hints, ann_hints):
def build_warnings(filter_stats, seqnames, ref_hints, ann_hints, ann=None):
"""Build all warning messages."""
warnings = []
if ann is not None and ann.was_gff and ann.only_in_gff:
examples = "\n".join(ann.only_in_gff[:5])
warnings.append(
f"Info: {len(ann.only_in_gff)} annotation "
"seqnames had no transcript records retained during GFF-to-GTF "
"conversion. This is expected for seqnames containing only "
"non-transcript GFF features.\n"
f"Examples:\n{examples}"
)
if filter_stats["excluded_unstranded_records"]:
warnings.extend([
UNSTRANDED_WARNING,
@ -413,15 +456,19 @@ def build_warnings(filter_stats, seqnames, ref_hints, ann_hints):
)
if seqnames["only_in_annotation"]:
warnings.extend([
"Warning: Some seqnames are present in the annotation but not the genome:",
*seqnames["only_in_annotation"][:5],
])
examples = "\n".join(seqnames["only_in_annotation"][:5])
warnings.append(
f"Warning: {len(seqnames['only_in_annotation'])} annotation "
"seqnames are not present in the reference.\n"
f"Examples:\n{examples}"
)
if seqnames["only_in_reference"]:
warnings.extend([
"Warning: Some seqnames are present in the genome but not the annotation:",
*seqnames["only_in_reference"][:5],
])
examples = "\n".join(seqnames["only_in_reference"][:5])
warnings.append(
f"Warning: {len(seqnames['only_in_reference'])} reference "
"seqnames are not present in the annotation.\n"
f"Examples:\n{examples}"
)
ref_builds = ref_hints["builds"]
ann_builds = ann_hints["builds"]
@ -466,7 +513,11 @@ def prepare_annotation_reference(annotation, reference, out_dir):
filter_stats = ann.filter_stranded()
# check seqname overlap between prepared annotation and reference
seqnames = validate_seqname_overlap(ann.seqnames, ref.seqnames)
seqnames = validate_seqname_overlap(
ann.seqnames,
ref.seqnames,
gff_seqnames=ann.only_in_gff,
)
if not seqnames["has_overlap"]:
raise ValueError(
@ -482,7 +533,7 @@ def prepare_annotation_reference(annotation, reference, out_dir):
# heuristically detect build and provider hints from filenames and headers
ref_hints = ref.detect_hints()
ann_hints = ann.detect_hints()
warnings = build_warnings(filter_stats, seqnames, ref_hints, ann_hints)
warnings = build_warnings(filter_stats, seqnames, ref_hints, ann_hints, ann)
if ann.normalised_records:
warnings.append(
"Warning: Replaced "
@ -500,6 +551,11 @@ def prepare_annotation_reference(annotation, reference, out_dir):
"normalised_unknown_strand_records": ann.normalised_records,
**filter_stats,
},
"gff_conversion": {
"input_count": ann.gff_convert_input_count,
"output_count": ann.gff_convert_output_count,
"pruned_count": ann.gff_convert_pruned_count,
},
"reference": {
"input": str(ref.input_path),
},
@ -508,6 +564,7 @@ def prepare_annotation_reference(annotation, reference, out_dir):
"seqname_overlap": seqnames["matches"],
"only_in_annotation": seqnames["only_in_annotation"],
"only_in_reference": seqnames["only_in_reference"],
"only_in_gff": ann.only_in_gff,
"reference_build_hints": ref_hints["builds"],
"annotation_build_hints": ann_hints["builds"],
"reference_provider_hints": ref_hints["providers"],

103
bin/workflow_glue/report.py
View File

@ -501,35 +501,20 @@ def _format_hint_values(hints):
def _create_warning_banner(message, level="warning"):
"""Create a styled warning banner."""
colors = {
"warning": "#fff3cd",
"danger": "#f8d7da",
"info": "#d1ecf1",
}
border_colors = {
"warning": "#ffc107",
"danger": "#dc3545",
"info": "#0dcaf0",
}
style = (
"padding: 15px; margin: 10px 0; "
"background-color: {}; "
"border-left: 4px solid {}; "
"border-radius: 4px;".format(
colors.get(level, colors["warning"]),
border_colors.get(level, border_colors["warning"]),
)
)
with div(style=style):
with strong():
raw(
alert_level = level if level in {"warning", "danger", "info"} else "warning"
icon = (
"⚠️"
if level == "warning"
if alert_level == "warning"
else "❌"
if level == "danger"
if alert_level == "danger"
else "ℹ️"
)
raw(message)
with div(
cls=f"alert alert-{alert_level}",
role="alert",
style="white-space: pre-line;",
):
raw(f"<strong>{icon}</strong> {escape(message)}")
def _heatmap_style():
@ -751,43 +736,91 @@ def main(args):
"combinations."
)
for warning in annotation_reference_summary.get("warnings", []):
_create_warning_banner(warning, level="warning")
_create_warning_banner(
warning,
level="info" if warning.startswith("Info:") else "warning",
)
seqname_rows = [
(
"Overlapping seqnames",
len(annotation_reference_summary.get("seqname_overlap", [])),
_format_count_value(
len(annotation_reference_summary.get("seqname_overlap", []))
),
),
(
"Seqnames only in annotation",
len(annotation_reference_summary.get("only_in_annotation", [])),
_format_count_value(
len(
annotation_reference_summary.get(
"only_in_annotation", []
)
)
),
),
(
"Seqnames only in reference",
len(annotation_reference_summary.get("only_in_reference", [])),
_format_count_value(
len(
annotation_reference_summary.get(
"only_in_reference", []
)
)
),
),
]
annotation_summary = annotation_reference_summary.get("annotation", {})
gff_conversion = annotation_reference_summary.get("gff_conversion", {})
if annotation_summary.get("was_gff"):
seqname_rows.extend([
(
"Annotation records retained",
annotation_summary.get("kept_records", "N/A"),
"Seqnames only in GFF (pruned during GTF conversion)",
_format_count_value(
len(annotation_reference_summary.get("only_in_gff", []))
),
),
])
annotation_rows = [
(
"Transcript annotation records retained",
_format_count_value(annotation_summary.get("kept_records")),
),
(
"Unstranded records excluded",
annotation_summary.get("excluded_unstranded_records", "N/A"),
_format_count_value(
annotation_summary.get("excluded_unstranded_records")
),
),
(
"Annotation attributes sanitised",
annotation_summary.get("sanitised_attribute_records", "N/A"),
"Transcript annotation attributes sanitised",
_format_count_value(
annotation_summary.get("sanitised_attribute_records")
),
])
),
]
if annotation_summary.get("was_gff"):
annotation_rows.append(
(
"Non-transcript records excluded "
"(pruned during GTF conversion)",
_format_count_value(gff_conversion.get("pruned_count")),
),
)
with div(cls="row"):
with div(cls="col-md-6"):
DataTable.from_pandas(
pd.DataFrame(seqname_rows, columns=["Check", "Value"]),
paging=False,
searchable=False,
use_index=False,
)
with div(cls="col-md-6"):
DataTable.from_pandas(
pd.DataFrame(annotation_rows, columns=["Check", "Value"]),
paging=False,
searchable=False,
use_index=False,
)
h4("Build and Provider Hints")
hint_rows = [

3
bin/workflow_glue/tests/common/test_report.py
View File

@ -247,7 +247,8 @@ def test_report_main_accepts_optional_file_sentinels(monkeypatch, tmp_path):
assert out_report.exists()
assert any("Overlapping seqnames" in table.to_string() for table in tables)
assert any(
"Annotation attributes sanitised" in table.to_string() for table in tables
"Transcript annotation attributes sanitised" in table.to_string()
for table in tables
)
assert any("GRCh38" in table.to_string() for table in tables)

2
bin/workflow_glue/tests/common/test_report_qc.py
View File

@ -168,7 +168,7 @@ def test_warning_banner_creation():
html = doc.render()
assert "Test warning message" in html
assert "background-color" in html
assert "alert alert-warning" in html
def test_bambu_qc_with_warnings():

64
bin/workflow_glue/tests/wf/test_prepare_annotation_reference.py
View File

@ -37,6 +37,13 @@ chr1\tsim\texon\t1\t5\t.\t+\t.\tParent=transcript1
chr1\tsim\texon\t6\t10\t.\t+\t.\tParent=transcript1
"""
GFF3_INPUT_WITH_REGION = """##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\t10\t.\t+\t.\tParent=transcript1
NT_167211.2\tRefSeq\tregion\t1\t176608\t.\t+\t.\tID=id1487055
"""
NCBI_REFERENCE = """>chr1 GRCh38 RefSeq primary assembly
AAAA
>chrExtra
@ -116,6 +123,14 @@ def test_prepare_all_formats_produce_valid_outputs(
# paths returned in summary must match
assert summary["annotation"]["prepared"] == str(out_dir / "annotation.gtf")
assert summary["reference"]["input"] == str(ref_path)
if summary["annotation"]["was_gff"]:
assert summary["gff_conversion"]["input_count"] is not None
assert summary["gff_conversion"]["output_count"] is not None
assert summary["gff_conversion"]["pruned_count"] is not None
else:
assert summary["gff_conversion"]["input_count"] is None
assert summary["gff_conversion"]["output_count"] is None
assert summary["gff_conversion"]["pruned_count"] is None
# outputs must be valid for downstream tools (Bambu, SQANTI)
reference = Path(summary["reference"]["input"])
@ -177,6 +192,55 @@ def test_gff3_conversion_via_gffread(tmp_path):
assert "\texon\t" in prepared_text
def test_gffread_omissions_are_reported_separately(tmp_path, monkeypatch):
"""GFF records and seqnames dropped by gffread are reported explicitly."""
reference = _write(
tmp_path / "reference.fa",
SIMPLE_REFERENCE + ">NT_167211.2\nAAAA\n",
)
annotation = _write(tmp_path / "annotation.gff3", GFF3_INPUT_WITH_REGION)
def fake_gffread_conversion(self, _gff_path, gtf_path):
Path(gtf_path).write_text(
'chr1\tsim\ttranscript\t1\t10\t.\t+\t.\t'
'gene_id "gene1"; transcript_id "transcript1";\n'
'chr1\tsim\texon\t1\t10\t.\t+\t.\t'
'gene_id "gene1"; transcript_id "transcript1";\n',
encoding="utf-8",
)
monkeypatch.setattr(
prepare_annotation_reference.Annotation,
"_run_gffread_conversion",
fake_gffread_conversion,
)
summary = prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
assert summary["gff_conversion"]["input_count"] == 4
assert summary["gff_conversion"]["output_count"] == 2
assert summary["gff_conversion"]["pruned_count"] == 2
assert summary["only_in_gff"] == ["NT_167211.2"]
assert "NT_167211.2" not in summary["seqnames"]["only_in_annotation"]
assert "NT_167211.2" not in summary["seqnames"]["only_in_reference"]
pruned_warning = next(
w
for w in summary["warnings"]
if "no transcript records retained during GFF-to-GTF conversion" in w
)
assert pruned_warning.startswith("Info:")
assert (
"expected for seqnames containing only non-transcript GFF features"
in pruned_warning
)
assert "Examples:\nNT_167211.2" in pruned_warning
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)