Decouple orf_class from ORF length across the custom/orf* modules

modules/nf-core/custom/orfcollapse/meta.yml

@@ -14,15 +14,17 @@ description: |
   Nat Biotechnol, doi:10.1038/s41587-022-01369-0; gencode-riboseqORFs
   collapse_cutoff 0.9), implemented here with MMseqs2 sequence-identity
   clustering rather than that tool's longest-shared-string / P-site-overlap
-  metric. Small ORFs (orf_class "smORF", i.e. aa_length <= 100) are clustered by
+  metric. Small ORFs (`aa_length` <= `--smorf-max-aa`, default 100) are clustered by
   amino-acid identity upstream and this module folds each multi-member cluster
   down to one representative.
 
-  Only smORF rows are collapsed; larger ORFs and transcript-anchored classes are
-  passed through untouched. Among the smORF members of a cluster the
-  representative is chosen by longest aa_length (ties broken by orf_id), so the
-  result does not depend on which sequence MMseqs2 labelled the cluster
-  representative. Catalogue row order is preserved; dropped members fold their
+  Only small ORFs are collapsed; larger ORFs are passed through untouched.
+  Eligibility is the catalogue's `is_smorf` flag, independent of `orf_class`, so
+  a short uORF and a short novel ORF are both candidates; `--smorf-max-aa`
+  re-derives the flag and aborts on disagreement. Among the members of a cluster
+  the representative is chosen by class specificity, then longest aa_length,
+  then orf_id, so the result does not depend on which sequence MMseqs2 labelled
+  the cluster representative. Catalogue row order is preserved; dropped members fold their
   `called_by_<caller>` / `score_<caller>` evidence, `n_samples` / `samples`
   recurrence and gene mappings into the survivor.
 keywords:

modules/nf-core/custom/orfcollapse/templates/orfcollapse.py

@@ -12,11 +12,11 @@
 two deliberate departures from that reference:
 
   - GENCODE collapses overlapping ORFs of any size within a shared locus;
-    this restricts collapsing to small ORFs (orf_class == "smORF", i.e.
-    aa_length <= 100) and clusters them locus-agnostically across the whole
+    this restricts collapsing to small ORFs (aa_length <= --smorf-max-aa,
+    default 100) and clusters them locus-agnostically across the whole
     catalogue, since the target case is one micropeptide recurring at several
-    non-overlapping loci. The smORF-only restriction is this pipeline's choice,
-    not a GENCODE property.
+    non-overlapping loci. The small-ORF-only restriction is this pipeline's
+    choice, not a GENCODE property.
   - similarity is MMseqs2 global sequence identity (--min-seq-id 0.9,
     mmseqs/easycluster upstream) rather than GENCODE's longest-shared-substring
     / P-site-overlap metric, so the 0.9 here approximates rather than
@@ -32,11 +32,14 @@
 high-confidence subset of the de-redundified catalogue and a folded
 micropeptide is judged on its combined cross-caller / cross-sample evidence.
 
-Only smORF rows are collapsed; larger ORFs and transcript-anchored classes pass
-through untouched, preserving the deterministic coordinate/transcript merge from
-upstream. Among the smORF members of a cluster the representative is chosen here
-(longest aa_length, ties broken by orf_id) so the result is independent of which
-sequence MMseqs2 labelled the cluster representative. Catalogue row order is
+Only small ORFs are collapsed; larger ORFs pass through untouched, preserving
+the deterministic coordinate/transcript merge from upstream. Eligibility is the
+catalogue's `is_smorf` flag, independent of `orf_class`, so a short uORF and a
+short novel ORF are both candidates; `--smorf-max-aa` re-derives the flag and
+aborts on disagreement. Among the members of a cluster the representative is
+chosen here (class specificity, then longest aa_length, then orf_id) so the
+result is independent of which sequence MMseqs2 labelled the cluster
+representative. Catalogue row order is
 preserved; dropped members fold their cross-caller / cross-sample evidence and
 gene mappings into the survivor.
 """
@@ -63,8 +66,12 @@
     "rpbp": "max",
     "price": "min",
 }
-CLASS_ORDER = ("canonical_cds", "uORF", "dORF", "novel_u", "smORF", "other")
-SMORF_CLASS = "smORF"
+CLASS_ORDER = ("canonical_cds", "uORF", "uoORF", "dORF", "doORF", "intORF", "novel_u", "other")
+
+# Survivor preference when a peptide cluster spans more than one class, most
+# specific first. Mirrors orfmerge's CLASS_SPECIFICITY so an annotated CDS is
+# never deleted by a longer novel ORF that happens to share its peptide.
+CLASS_SPECIFICITY = ("canonical_cds", "uoORF", "uORF", "doORF", "dORF", "intORF", "novel_u", "other")
 
 
 def read_fasta(path):
@@ -113,15 +120,26 @@ def best_score(values, direction):
 
 
 def merge_members(members):
-    """Fold smORF rows sharing an AA cluster into one representative row dict."""
-    rep = sorted(members, key=lambda r: (-int(r.get("aa_length") or 0), r["orf_id"]))[0]
+    """Fold small-ORF rows sharing an AA cluster into one representative row dict."""
+    rank = {c: i for i, c in enumerate(CLASS_SPECIFICITY)}
+    rep = sorted(
+        members,
+        key=lambda r: (
+            rank.get(r.get("orf_class", "other"), len(rank)),
+            -int(r.get("aa_length") or 0),
+            r["orf_id"],
+        ),
+    )[0]
     out = dict(rep)
     for c in CALLERS:
         out[f"called_by_{c}"] = "1" if any(r.get(f"called_by_{c}") == "1" for r in members) else "0"
         out[f"score_{c}"] = best_score([r.get(f"score_{c}", "") for r in members], SCORE_DIRECTIONS[c])
     samples = sorted({s for r in members for s in (r.get("samples") or "").split(",") if s})
     out["n_samples"] = str(len(samples))
     out["samples"] = ",".join(samples)
+    if "orf_type_native" in out:
+        natives = sorted({t for r in members for t in (r.get("orf_type_native") or "").split(",") if t})
+        out["orf_type_native"] = ",".join(natives)
     return rep["orf_id"], out
 
 
@@ -139,14 +157,36 @@ def main():
         default=1,
         help="Minimum distinct samples for an ORF to enter the consensus view (default: 1)",
     )
+    parser.add_argument(
+        "--smorf-max-aa",
+        type=int,
+        default=100,
+        help="Maximum aa_length eligible for peptide-level collapse (default: 100)",
+    )
     args = parser.parse_args(shlex.split("${args}"))
+    if args.smorf_max_aa < 1:
+        sys.exit(f"orfcollapse: --smorf-max-aa must be >= 1, got {args.smorf_max_aa}")
 
     prefix = "${prefix}"
 
     catalogue = pd.read_csv("${catalogue_tsv}", sep="\\t", comment="#", dtype=str, keep_default_na=False)
     header = list(catalogue.columns)
     rows = catalogue.to_dict("records")
 
+    # The collapse scope is derived from these columns, so a silent rename
+    # upstream must abort rather than quietly collapse nothing.
+    missing = [c for c in ("orf_class", "aa_length", "is_smorf") if c not in header]
+    if missing:
+        sys.exit(f"orfcollapse: catalogue is missing required column(s) {missing}")
+
+    unknown = sorted(set(catalogue["orf_class"]) - set(CLASS_ORDER))
+    if unknown:
+        sys.exit(f"orfcollapse: unknown orf_class value(s) {unknown}; update CLASS_ORDER")
+
+    invalid_smorf = sorted(set(catalogue["is_smorf"]) - {"0", "1"})
+    if invalid_smorf:
+        sys.exit(f"orfcollapse: invalid is_smorf value(s) {invalid_smorf}; expected 0 or 1")
+
     bed_index = {}
     with open("${bed12}") as fh:
         for line in fh:
@@ -157,9 +197,26 @@ def main():
     aa = read_fasta("${aa_fasta}")
     cluster_of = read_clusters("${cluster_tsv}")
 
+    # Eligibility is the catalogue's is_smorf; --smorf-max-aa is kept only to
+    # re-derive it, so a threshold mismatch between processes fails loudly.
+    def is_small(row):
+        try:
+            aa = int(row.get("aa_length") or 0)
+        except ValueError:
+            return False
+        return 0 < aa <= args.smorf_max_aa
+
+    diverged = [r["orf_id"] for r in rows if (r.get("is_smorf") == "1") != is_small(r)]
+    if diverged:
+        sys.exit(
+            f"orfcollapse: is_smorf disagrees with --smorf-max-aa={args.smorf_max_aa} on "
+            f"{len(diverged)} row(s) (first: {diverged[:3]}); pass the same --smorf-max-aa "
+            "to CUSTOM_ORFNORMALISE and CUSTOM_ORFCOLLAPSE"
+        )
+
     clusters = defaultdict(list)
     for r in rows:
-        if r.get("orf_class") == SMORF_CLASS:
+        if r.get("is_smorf") == "1":
             clusters[cluster_of.get(r["orf_id"], r["orf_id"])].append(r)
 
     remap, merged_rows, dropped = {}, {}, set()

modules/nf-core/custom/orfcollapse/tests/main.nf.test

@@ -9,6 +9,8 @@ nextflow_process {
     tag "custom"
     tag "custom/orfcollapse"
 
+    config "./nextflow.config"
+
     test("collapse identical small ORFs across loci") {
 
         when {
@@ -38,6 +40,9 @@ nextflow_process {
                 { assert process.success },
                 { assert snapshot(sanitizeOutput(process.out)).match() },
                 // The two identical smORFs on different strands/loci collapse to one row.
+                // Four input rows fold to three, and the assertions below pin
+                // which row went and that its evidence moved to the survivor,
+                // so a silently empty eligible set cannot pass unnoticed.
                 { assert rows.size() == 3 },
                 { assert byId.containsKey('orf_00000001') },
                 { assert !byId.containsKey('orf_00000002') },
@@ -86,4 +91,34 @@ nextflow_process {
         }
     }
 
+
+    test("a --smorf-max-aa that disagrees with is_smorf aborts") {
+
+        when {
+            params {
+                module_args = "--smorf-max-aa 10"
+            }
+            process {
+                """
+                def dir = params.modules_testdata_base_path + "genomics/homo_sapiens/riboseq_expression/orf_catalogue"
+                input[0] = [
+                    [ id: 'cohort' ],
+                    file("\${dir}/cohort.catalogue.bed12",   checkIfExists: true),
+                    file("\${dir}/cohort.catalogue.tsv",     checkIfExists: true),
+                    file("\${dir}/cohort.orf_to_gene.tsv",   checkIfExists: true),
+                    file("\${dir}/cohort.catalogue.aa.fasta",checkIfExists: true),
+                    file("\${dir}/cohort_cluster.tsv",       checkIfExists: true)
+                ]
+                """
+            }
+        }
+
+        then {
+            assertAll(
+                { assert process.failed },
+                { assert process.errorReport.contains('is_smorf disagrees') }
+            )
+        }
+    }
+
 }

modules/nf-core/custom/orfcollapse/tests/main.nf.test.snap

@@ -23,7 +23,7 @@
                         {
                             "id": "cohort"
                         },
-                        "cohort.catalogue.tsv:md5,2ae2eb7006ab613d0d23ea5eaa959615"
+                        "cohort.catalogue.tsv:md5,a9007ad196e22f8817576368c9036bbf"
                     ]
                 ],
                 "consensus_bed12": [
@@ -47,15 +47,15 @@
                         {
                             "id": "cohort"
                         },
-                        "cohort.catalogue.consensus.tsv:md5,2ae2eb7006ab613d0d23ea5eaa959615"
+                        "cohort.catalogue.consensus.tsv:md5,a9007ad196e22f8817576368c9036bbf"
                     ]
                 ],
                 "multiqc": [
                     [
                         {
                             "id": "cohort"
                         },
-                        "cohort.catalogue.mqc.tsv:md5,0541bc7510de3bc238145a4f913a9be9"
+                        "cohort.catalogue.mqc.tsv:md5,29287249a0db954902f07adaa72b6bb4"
                     ]
                 ],
                 "orf_to_gene_tsv": [
@@ -71,11 +71,11 @@
                 ]
             }
         ],
+        "timestamp": "2026-08-03T10:46:50.930367",
         "meta": {
-            "nf-test": "0.9.3",
-            "nextflow": "25.04.8"
-        },
-        "timestamp": "2026-06-25T20:22:39.073880596"
+            "nf-test": "0.9.5",
+            "nextflow": "25.10.4"
+        }
     },
     "collapse identical small ORFs across loci - stub": {
         "content": [
@@ -149,10 +149,10 @@
                 ]
             }
         ],
+        "timestamp": "2026-08-03T10:46:57.809342",
         "meta": {
-            "nf-test": "0.9.3",
-            "nextflow": "25.04.8"
-        },
-        "timestamp": "2026-06-25T20:22:43.127779465"
+            "nf-test": "0.9.5",
+            "nextflow": "25.10.4"
+        }
     }
 }

modules/nf-core/custom/orfcollapse/tests/nextflow.config

@@ -0,0 +1,7 @@
+params.module_args = ""
+
+process {
+    withName: "CUSTOM_ORFCOLLAPSE" {
+        ext.args = params.module_args
+    }
+}

modules/nf-core/custom/orfmerge/meta.yml

@@ -6,19 +6,29 @@ description: |
   (typically) `bedtools/getfasta` + `seqkit/translate` downstream to obtain
   the AA FASTA.
 
-  Strategy is class-aware (operating on the harmonised `orf_class` written
-  by `custom/orfnormalise`):
+  Rows sharing an identical exon structure are grouped first and clustered as a
+  single proxy chosen by class specificity, then the members are restored. Class
+  selects the strategy, so without this an ORF that two callers class differently
+  would be split across strategies and could never merge.
 
-    - canonical_cds:        collapse by (transcript_id, strand). One
-                            canonical CDS per transcript by definition.
-    - uORF, dORF, other:    collapse by (transcript_id, strand, start,
-                            end). A single transcript can host multiple
+  Proxies are then partitioned by clustering strategy. The harmonised `orf_class`
+  written by `custom/orfnormalise` selects the strategy but is never part of a
+  grouping key, because callers disagree on class for the same ORF (Ribo-TISH
+  reports `5'UTR` for both uORFs and CDS-overlapping uORFs) and keying on it
+  would emit one catalogue row per disagreeing caller:
+
+    - canonical_cds:        grouped by (transcript_id, strand), then
+                            reciprocal-overlap clustered within the
+                            transcript so a short truncated variant is
+                            not folded into the full-length CDS.
+    - uORF, uoORF, dORF,    collapse by (transcript_id, strand, start,
+      doORF, intORF, other: end). A single transcript can host multiple
                             distinct uORFs / dORFs / internal ORFs, so
                             keying on the outer span keeps them in
                             separate clusters while still merging
                             cross-caller calls that agree on coordinates.
-    - novel_u, smORF:       greedy reciprocal-overlap clustering on the
-                            outer genomic span at `--reciprocal-overlap`
+    - novel_u:              greedy reciprocal-overlap clustering on summed
+                            exon-block intersection at `--reciprocal-overlap`
                             (default 0.8). Catches fuzzy cross-caller
                             matches and exact-coordinate collapses in
                             one pass. Order-dependent at the boundary:
@@ -112,7 +122,9 @@ output:
             Per-ORF table with `called_by_<caller>` and `score_<caller>`
             columns for each supported caller, `n_samples` / `samples`
             cross-sample recurrence columns, plus orf_class, aa_length,
-            and host gene_id / transcript_id.
+            host gene_id / transcript_id, `orf_type_native` listing every
+            caller-native ORF-type label in the cluster, and the
+            `is_smorf` length flag.
           pattern: "*.tsv"
           ontologies:
             - edam: http://edamontology.org/format_3475 # TSV