Merge branch 'v202' into 'dev'

wf-transcriptomes v2.0.2

See merge request epi2melabs/workflows/wf-transcriptomes!348

.gitlab-ci.yml

@@ -191,6 +191,13 @@ docker-run:
                 --ref_annotation ${CI_PROJECT_NAME}/data/gencode.v22.annotation.chr20.gtf \
                 --sample_sheet ${CI_PROJECT_NAME}/data/mods_rna_subset/sample_sheet.csv \
                 --igv"
+            AFTER_NEXTFLOW_CMD: >
+                for sample in sample01 sample02; do
+                  for suffix in mods.bedmethyl.gz mods.inosine.bw mods.m5C.bw mods.m6A.bw mods.pseU.bw mods.summary.tsv; do
+                    test -f ${CI_PROJECT_NAME}/samples/$${sample}/mods/$${sample}.$${suffix};
+                  done;
+                done;
+                
         - if: $MATRIX_NAME == "mouse_splice_error"
           variables:
               NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/mouse_splice_fail.tar.gz -O ${CI_PROJECT_NAME}/data/mouse_splice_fail.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/mouse_splice_fail.tar.gz  -C ${CI_PROJECT_NAME}/data/" 

CHANGELOG.md

@@ -5,11 +5,23 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 
-## [Unreleased]
+## [v2.0.2]
+
+This patch release of `wf-transcriptomes` fixes the `minimap2` invocation for cDNA read mapping, ensuring that both strands are searched for canonical splice sites.
+This release also distinguishes annotation seqnames with no transcript records after GFF-to-GTF conversion from seqnames that are absent from the annotation.
+Users analysing cDNA data should adopt this release.
 
 ### 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.
+- Dataframe merge error during PCA plotting caused by all numeric aliases in the sample sheet.
+
+### 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".
+
+### Changed
+- 2Ome* mod code labels are used in place of CHEBI numbers in output file names and reports for 2'-O-methylation modifications.
+- Removed mention of `analysis_group` and `type` sample sheet columns from `README.md` as they are not relevant for this workflow.
+
 
 ## [v2.0.1]
 

bin/workflow_glue/prepare_annotation_reference.py

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

bin/workflow_glue/report.py

@@ -63,6 +63,22 @@ classification_categories = {
     "Intergenic": "Query isoform lies in an intergenic region.",
 }
 
+transcript_meta_dtypes = {
+    "TXNAME": "string",
+    "GENEID": "string",
+    "NDR": "Float64",
+    "novelGene": "boolean",
+    "novelTranscript": "boolean",
+    "txClassDescription": "string",
+    "readCount": "Int64",
+    "relReadCount": "Float64",
+    "relSubsetCount": "Float64",
+    "txid": "string",
+    "eqClassById": "string",
+    "gene_name": "string",
+    "transcript_name": "string"
+}
+
 
 def get_bokeh_widgets_js():
     """Return the inline Bokeh widgets JavaScript bundle."""
@@ -114,20 +130,15 @@ def _format_ratio_value(value):
 
 def _sorted_transcript_abundance_table(tx_counts_file, sample_aliases):
     """Get transcript abundance rows by descending total abundance of sample columns."""
-    tx_counts = _read_table(tx_counts_file)
+    dtype = transcript_meta_dtypes | {sample: "Float64" for sample in sample_aliases}
+    tx_counts = _read_table(tx_counts_file, dtype=dtype)
+
     if tx_counts is None or tx_counts.empty:
         return tx_counts
-
-    count_columns = [
-        column for column in tx_counts.columns if column in set(sample_aliases)
-    ]
-    if not count_columns:
-        return tx_counts
-
-    numeric_counts = tx_counts[count_columns].apply(pd.to_numeric, errors="coerce")
-    abundance = numeric_counts.sum(axis=1, min_count=1)
-    order = abundance.fillna(float("-inf")).sort_values(ascending=False).index
-    return tx_counts.loc[order].reset_index(drop=True)
+    # Single sample data will not have 'readCount' column,
+    # so sort by the first sample alias counts instead.
+    sort_col = 'readCount' if 'readCount' in tx_counts.columns else sample_aliases[0]
+    return tx_counts.sort_values(sort_col, ascending=False)
 
 
 def _format_classification_label(name):
@@ -137,7 +148,10 @@ def _format_classification_label(name):
 
 def _transcriptome_summary(transcriptome_dir):
     """Return transcriptome metrics and transcript class counts DataFrames."""
-    tx_meta = _read_table(Path(transcriptome_dir) / "transcript_metadata.tsv")
+    tx_meta = _read_table(
+        Path(transcriptome_dir) / "transcript_metadata.tsv",
+        dtype=transcript_meta_dtypes
+    )
     if tx_meta is None:
         return None, None
 
@@ -168,7 +182,10 @@ def _sample_summaries(samples_dir):
     for sample_dir in sorted(samples_path.iterdir()):
         if not sample_dir.is_dir():
             continue
-        tx_meta = _read_table(sample_dir / "transcript_metadata.tsv")
+        tx_meta = _read_table(
+            sample_dir / "transcript_metadata.tsv",
+            dtype=transcript_meta_dtypes
+        )
         if tx_meta is None:
             continue
         summaries[sample_dir.name] = pd.DataFrame(
@@ -196,7 +213,13 @@ def _sqanti_table(sqanti_dir):
         if "classification_summary.tsv" in files:
             summaries.append(Path(root) / "classification_summary.tsv")
     for summary in sorted(summaries):
-        table = _read_table(summary)
+        table = _read_table(
+            summary,
+            dtype={
+                'structural_category': 'string',
+                'count': 'Int64'
+            }
+        )
         if table is None or table.empty:
             continue
 
@@ -235,7 +258,17 @@ def _sample_mod_summaries(summary_dir):
     for summary_file in sorted(summaries_path.glob("*.mods.summary.tsv")):
         if not summary_file.is_file():
             continue
-        summary = _read_table(summary_file)
+        summary = _read_table(
+            summary_file,
+            dtype={
+                "sample": "string",
+                "full_mod_code": "string",
+                "mod_label": "string",
+                "valid_coverage": "Int64",
+                "modified_calls": "Int64",
+                "modification_percent": "Float64"
+            },
+        )
         if summary is None or summary.empty:
             continue
         if "sample" not in summary.columns:
@@ -395,7 +428,7 @@ def _contrast_results(de_dir, filename):
         # Enforce str dtype in case of all Nan values.
         table = _read_table(
             contrast_dir / filename,
-            dtype={'gene_name': str, 'transcript_name': str}
+            dtype={'gene_name': 'string', 'transcript_name': 'string'}
         )
         if table is None or table.empty:
             continue
@@ -466,15 +499,32 @@ def _load_de_qc(de_dir):
     return None
 
 
-def _load_cpm_tables(cohort_dir):
+def _load_cpm_tables(cohort_dir, sample_aliases):
     """Load cohort-level gene and transcript CPM tables."""
     cohort_dir = Path(cohort_dir)
-    gene_cpm = _read_table(cohort_dir / "gene_cpm.tsv")
+    gene_columns = ["GENEID", *sample_aliases]
+    transcript_columns = ["TXNAME", *sample_aliases]
+
+    gene_cpm = _read_table(
+        cohort_dir / "gene_cpm.tsv",
+        dtype={
+            "GENEID": "string",
+            **{sample: "Float64" for sample in sample_aliases},
+        },
+        usecols=gene_columns,
+    )
 
     if gene_cpm is None or gene_cpm.empty:
         gene_cpm = None
 
-    transcript_cpm = _read_table(cohort_dir / "transcript_cpm.tsv")
+    transcript_cpm = _read_table(
+        cohort_dir / "transcript_cpm.tsv",
+        dtype={
+            "TXNAME": "string",
+            **{sample: "Float64" for sample in sample_aliases},
+        },
+        usecols=transcript_columns,
+    )
     if transcript_cpm is None or transcript_cpm.empty:
         transcript_cpm = None
 
@@ -489,7 +539,7 @@ def _load_cohort_samples(cohort_dir):
     sample_file = Path(cohort_dir) / "samples.csv"
     if not sample_file.exists():
         return None
-    return pd.read_csv(sample_file)
+    return pd.read_csv(sample_file, dtype='string')
 
 
 def _format_hint_values(hints):
@@ -501,35 +551,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"]),
-        )
+    alert_level = level if level in {"warning", "danger", "info"} else "warning"
+    icon = (
+        "⚠️"
+        if alert_level == "warning"
+        else "❌"
+        if alert_level == "danger"
+        else "ℹ️"
     )
-    with div(style=style):
-        with strong():
-            raw(
-                "⚠️ "
-                if level == "warning"
-                else "❌ "
-                if 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():
@@ -662,6 +697,7 @@ def main(args):
 
     with open(args.metadata, "r") as handle:
         metadata = json.load(handle)
+    sample_aliases = [item["alias"] for item in metadata]
 
     if args.stats:
         with report.add_section("Read summary", "Reads"):
@@ -751,43 +787,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", {})
-            seqname_rows.extend([
+            gff_conversion = annotation_reference_summary.get("gff_conversion", {})
+            if annotation_summary.get("was_gff"):
+                seqname_rows.extend([
+                    (
+                        "Seqnames only in GFF (pruned during GTF conversion)",
+                        _format_count_value(
+                            len(annotation_reference_summary.get("only_in_gff", []))
+                        ),
+                    ),
+                ])
+            annotation_rows = [
                 (
-                    "Annotation records retained",
-                    annotation_summary.get("kept_records", "N/A"),
+                    "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")
+                    ),
                 ),
-            ])
-            DataTable.from_pandas(
-                pd.DataFrame(seqname_rows, columns=["Check", "Value"]),
-                paging=False,
-                searchable=False,
-                use_index=False,
-            )
+            ]
+            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 = [
@@ -1005,7 +1089,7 @@ def main(args):
 
         tx_counts = _sorted_transcript_abundance_table(
             Path(bambu_dir) / "transcript_counts.tsv",
-            [item["alias"] for item in metadata]
+            sample_aliases
         )
         if tx_counts is not None and not tx_counts.empty:
             p(
@@ -1349,7 +1433,7 @@ def main(args):
 
         if de_qc:
             condition_column = de_qc.get("condition_column")
-            cohort_cpm = _load_cpm_tables(args.cohort_dir)
+            cohort_cpm = _load_cpm_tables(args.cohort_dir, sample_aliases)
             cohort_samples = _load_cohort_samples(args.cohort_dir)
 
         with report.add_section("Differential gene expression", "DGE"):

bin/workflow_glue/tests/common/test_report.py

@@ -45,9 +45,21 @@ def _write(path, text):
 
 def _build_report_args(tmp_path, de_qc=None):
     """Create minimal report inputs, optionally including DE QC JSON."""
+    # When DE/DTU QC inputs are present, the report also loads
+    # cohort CPM tables and samples.csv, so metadata.json must
+    # match the sample aliases in those files.
+    metadata_rows = [{"alias": "sampleA", "has_stats": False}]
+    if de_qc is not None and "contrasts" in de_qc:
+        metadata_rows = [
+            {"alias": f"sample_control_{i}", "has_stats": False}
+            for i in range(3)
+        ] + [
+            {"alias": f"sample_treated_{i}", "has_stats": False}
+            for i in range(3)
+        ]
     metadata = _write(
         tmp_path / "metadata.json",
-        json.dumps([{"alias": "sampleA", "has_stats": False}]),
+        json.dumps(metadata_rows),
     )
     params = _write(tmp_path / "params.json", "{}")
     versions = tmp_path / "versions"
@@ -247,7 +259,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)
 
@@ -797,15 +810,15 @@ def test_round_de_table_uses_scientific_notation_when_fixed_decimal_would_zero()
     assert rounded.loc[0, "padj"] == "4.900e-04"
 
 
-def test_sorted_transcript_abundance_table_uses_sample_alias_columns(tmp_path):
-    """Transcript abundance sorting should use only real sample alias columns."""
+def test_sorted_transcript_abundance_table_sorting(tmp_path):
+    """Transcript abundance sorting should sort by sample cohort readCount."""
     tx_counts_file = _write(
         tmp_path / "transcript_counts.tsv",
         (
-            "TXNAME\tsampleA\tsampleB\tgene_length\tannotation_score\n"
-            "tx_low\t1\t1\t10000\t999\n"
-            "tx_high\t5\t5\t10\t1\n"
-            "tx_mid\t2\t2\t5000\t500\n"
+            "TXNAME\treadCount\tsampleA\tsampleB\tgene_length\tannotation_score\n"
+            "tx_low\t2\t1\t1\t10000\t999\n"
+            "tx_high\t8\t5\t5\t10\t1\n"
+            "tx_mid\t4\t2\t2\t5000\t500\n"
         ),
     )
 
@@ -814,3 +827,23 @@ def test_sorted_transcript_abundance_table_uses_sample_alias_columns(tmp_path):
     )
 
     assert sorted_table["TXNAME"].tolist() == ["tx_high", "tx_mid", "tx_low"]
+
+
+def test_read_table_preserves_na_with_nullable_float_dtype(tmp_path):
+    """_read_table should preserve NA values with pandas nullable Float64."""
+    table_file = _write(
+        tmp_path / "nullable_float.tsv",
+        (
+            "name\tvalue\n"
+            "row1\t1.5\n"
+            "row2\tNA\n"
+        ),
+    )
+
+    table = report._read_table(
+        table_file, dtype={"value": "Float64"}, index_col="name"
+    )
+
+    assert str(table["value"].dtype) == "Float64"
+    assert table.loc["row1", "value"] == 1.5
+    assert pd.isna(table.loc["row2", "value"])

bin/workflow_glue/tests/common/test_report_qc.py

@@ -129,7 +129,7 @@ def test_load_cpm_tables_exists(tmp_path):
         "tx1\t3.0\t4.0\n"
     )
 
-    result = _load_cpm_tables(cohort_dir)
+    result = _load_cpm_tables(cohort_dir, sample_aliases=["sample1", "sample2"])
     assert result["gene"] is not None
     assert result["transcript"] is not None
     assert list(result["gene"]["GENEID"]) == ["gene1"]
@@ -143,7 +143,7 @@ def test_load_cpm_tables_missing(tmp_path):
     cohort_dir = tmp_path / "cohort"
     cohort_dir.mkdir()
 
-    result = _load_cpm_tables(cohort_dir)
+    result = _load_cpm_tables(cohort_dir, sample_aliases=[])
     assert result["gene"] is None
     assert result["transcript"] is None
 
@@ -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():

bin/workflow_glue/tests/wf/test_prepare_annotation_reference.py

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

data/mod_code_labels.tsv

@@ -3,3 +3,7 @@ A:17596	inosine
 A:a	m6A
 C:m	m5C
 T:17802	pseU
+C:19228	2OmeC
+A:69426	2OmeA
+G:19229	2OmeG
+T:19227	2OmeU

nextflow.config

@@ -67,7 +67,7 @@ manifest {
     description     = 'Long-read transcript discovery, quantification, differential expression, QC and mod counting.'
     mainScript      = 'main.nf'
     nextflowVersion = '>=23.04.2'
-    version         = 'v2.0.1'
+    version         = 'v2.0.2'
 }
 
 process {