-
Notifications
You must be signed in to change notification settings - Fork 1.1k
Decouple orf_class from ORF length across the custom/orf* modules #12498
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Changes from 6 commits
baab83c
ef4245e
9875b43
9f79e35
b8a0b6e
ac4d9be
80a89eb
23abd3f
59f659e
90adc39
d121556
05f9d4d
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change | ||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
@@ -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,10 +32,11 @@ | |||||||||||||||||||||||||||||||||||||||||
| 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 | ||||||||||||||||||||||||||||||||||||||||||
| Only small ORFs are collapsed; larger ORFs pass through untouched, preserving | ||||||||||||||||||||||||||||||||||||||||||
| the deterministic coordinate/transcript merge from upstream. Eligibility is a | ||||||||||||||||||||||||||||||||||||||||||
| length test, independent of `orf_class`, so a short uORF and a short novel ORF | ||||||||||||||||||||||||||||||||||||||||||
| are both candidates. Among the 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 | ||||||||||||||||||||||||||||||||||||||||||
| preserved; dropped members fold their cross-caller / cross-sample evidence and | ||||||||||||||||||||||||||||||||||||||||||
| gene mappings into the survivor. | ||||||||||||||||||||||||||||||||||||||||||
|
|
@@ -63,8 +64,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 +118,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 +155,32 @@ 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}") | ||||||||||||||||||||||||||||||||||||||||||
|
Comment on lines
+160
to
+168
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [Medium] Remove the duplicate smORF threshold
Suggested change
Comment added by Codex
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Kept the flag, but made the disagreement impossible rather than removing the cross-check.
The reason for keeping it: deleting the argument removes the disagreement but also removes the check, and the module stays usable standalone outside the pipeline. Happy to drop it entirely if you'd rather orfcollapse take the catalogue's word unconditionally. |
||||||||||||||||||||||||||||||||||||||||||
|
|
||||||||||||||||||||||||||||||||||||||||||
| 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 two columns, so a silent rename | ||||||||||||||||||||||||||||||||||||||||||
| # upstream must abort rather than quietly collapse nothing. | ||||||||||||||||||||||||||||||||||||||||||
| missing = [c for c in ("orf_class", "aa_length") 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") | ||||||||||||||||||||||||||||||||||||||||||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [Medium] Require and validate Collapse now depends on the flag from the merged catalogue, so it should fail clearly if that column is missing or contains anything except
Suggested change
Comment added by Codex
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Done — |
||||||||||||||||||||||||||||||||||||||||||
|
|
||||||||||||||||||||||||||||||||||||||||||
| bed_index = {} | ||||||||||||||||||||||||||||||||||||||||||
| with open("${bed12}") as fh: | ||||||||||||||||||||||||||||||||||||||||||
| for line in fh: | ||||||||||||||||||||||||||||||||||||||||||
|
|
@@ -157,9 +191,18 @@ def main(): | |||||||||||||||||||||||||||||||||||||||||
| aa = read_fasta("${aa_fasta}") | ||||||||||||||||||||||||||||||||||||||||||
| cluster_of = read_clusters("${cluster_tsv}") | ||||||||||||||||||||||||||||||||||||||||||
|
|
||||||||||||||||||||||||||||||||||||||||||
| # Eligibility is derived here from aa_length rather than read from a | ||||||||||||||||||||||||||||||||||||||||||
| # propagated flag, so the scope cannot silently drift from --smorf-max-aa. | ||||||||||||||||||||||||||||||||||||||||||
| def is_small(row): | ||||||||||||||||||||||||||||||||||||||||||
| try: | ||||||||||||||||||||||||||||||||||||||||||
| aa = int(row.get("aa_length") or 0) | ||||||||||||||||||||||||||||||||||||||||||
| except ValueError: | ||||||||||||||||||||||||||||||||||||||||||
| return False | ||||||||||||||||||||||||||||||||||||||||||
| return 0 < aa <= args.smorf_max_aa | ||||||||||||||||||||||||||||||||||||||||||
|
|
||||||||||||||||||||||||||||||||||||||||||
| clusters = defaultdict(list) | ||||||||||||||||||||||||||||||||||||||||||
| for r in rows: | ||||||||||||||||||||||||||||||||||||||||||
| if r.get("orf_class") == SMORF_CLASS: | ||||||||||||||||||||||||||||||||||||||||||
| if is_small(r): | ||||||||||||||||||||||||||||||||||||||||||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [Medium] Use the flag already in the catalogue Use the decision made by
Suggested change
Comment added by Codex
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Done — |
||||||||||||||||||||||||||||||||||||||||||
| clusters[cluster_of.get(r["orf_id"], r["orf_id"])].append(r) | ||||||||||||||||||||||||||||||||||||||||||
|
|
||||||||||||||||||||||||||||||||||||||||||
| remap, merged_rows, dropped = {}, {}, set() | ||||||||||||||||||||||||||||||||||||||||||
|
|
||||||||||||||||||||||||||||||||||||||||||
Uh oh!
There was an error while loading. Please reload this page.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[Low] Update the docs to match the code
The docs say that the longest ORF wins, but the code uses class specificity first. They should also name
is_smorfas the eligibility source.Comment added by Codex
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done —
orfcollapse/meta.ymlnow reads "the representative is chosen by class specificity, then longest aa_length, then orf_id", and names the catalogue'sis_smorfflag as the eligibility source.