From d4e11d6dda62329fa481ab90685d152191c09a9e Mon Sep 17 00:00:00 2001 From: aliyz Date: Wed, 30 Sep 2026 17:05:04 +0200 Subject: [PATCH 1/5] feat: protein trace events --- FastOMA.nf | 33 ++++++ FastOMA/_infer_subhog.py | 27 +++++ FastOMA/_protein_trace.py | 82 +++++++++++++++ FastOMA/_utils_roothog.py | 35 +++++- FastOMA/protein_trace.py | 216 ++++++++++++++++++++++++++++++++++++++ nextflow.config | 1 + nextflow_schema.json | 5 + pyproject.toml | 1 + 8 files changed, 399 insertions(+), 1 deletion(-) create mode 100644 FastOMA/_protein_trace.py create mode 100644 FastOMA/protein_trace.py diff --git a/FastOMA.nf b/FastOMA.nf index e3a27d2..5c3b6b3 100644 --- a/FastOMA.nf +++ b/FastOMA.nf @@ -146,6 +146,7 @@ process infer_roothogs{ path "omamer_rhogs" path "gene_id_dic_xml.pickle" path "selected_isoforms" , optional: true + path "protein_events_*.events", optional: true, emit: protein_events script: """ @@ -209,6 +210,7 @@ process hog_big{ path "*.nwk" , optional: true // gene trees if write True path "*.tsv", optional: true path "*.tsv.gz", optional: true + path "protein_events_*.events", optional: true, emit: protein_events script: """ @@ -247,6 +249,7 @@ process hog_rest{ path "*.nwk" , optional: true // gene trees if write True path "*.tsv", optional: true path "*.tsv.gz", optional: true + path "protein_events_*.events", optional: true, emit: protein_events script: """ fastoma-infer-subhogs --input-rhog-folder ${rhogsrest} \ @@ -299,6 +302,28 @@ process collect_subhogs{ """ } +process protein_trace { + label "process_single" + publishDir params.output_folder, mode: 'copy' + + input: + path proteome_folder + path events, stageAs: "events/*" + path roothogs_tsv + + output: + path "protein_fates.tsv" + path "protein_events.tsv.gz" + path "protein_trace_summary.txt" + + script: + """ + python -m FastOMA.protein_trace --proteomes ${proteome_folder} \ + --events events/ \ + --roothog-tsv ${roothogs_tsv} + """ +} + process extract_pairwise_ortholog_relations { label "process_medium" publishDir params.output_folder, mode: 'copy' @@ -492,6 +517,14 @@ workflow { channel.empty().concat(pickle_big_rhog, pickle_rest_rhog).set{ all_rhog_pickle } (orthoxml_file, OrthologousGroupsFasta, OrthologousGroups_tsv, rootHOGs_tsv) = collect_subhogs(all_rhog_pickle.collect(), gene_id_dic_xml, omamer_rhogs, species_tree_checked, params.fasta_header_id_transformer) + if (params.protein_trace) { + channel.empty() + .concat(infer_roothogs.out.protein_events, hog_big.out.protein_events, hog_rest.out.protein_events) + .flatten() + .collect() + .set { all_protein_events } + protein_trace(proteome_folder, all_protein_events, rootHOGs_tsv) + } c = hogmap.count().branch{ n-> TRUE: (n<=25 || params.force_pairwise_ortholog_generation) FALSE: n>25 diff --git a/FastOMA/_infer_subhog.py b/FastOMA/_infer_subhog.py index ad14bd6..f49855e 100644 --- a/FastOMA/_infer_subhog.py +++ b/FastOMA/_infer_subhog.py @@ -26,6 +26,7 @@ from . import _wrappers, logger from . import _utils_subhog from . import _utils_frag_SO_detection +from . import _protein_trace as ptrace from ._hog_class import HOG, Representative, split_hog, attach_scores, _member_species, _species_name_index, ScoreFlags from ._utils_subhog import MSAFilter, MSAFilterElbow, MSAFilterTrimAL @@ -104,8 +105,13 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder, pi tot_genes, placed_genes = 0, 0 hogs_rhogs_xml = [] + members_final = set() for hog_i in hogs_a_rhog: tot_genes += len(hog_i) + members_final |= {part for m in hog_i.get_members() for part in m.split("_|_")} # merged fragments count per protein + final_event = "in_final_hog" if len(hog_i) >= inferhog_min_hog_size_xml else "final_singleton_not_written" + for member in hog_i.get_members(): + ptrace.member_event(member, "infer_subhogs", final_event, rhogid, f"hog={hog_i.hogid} size={len(hog_i)}") if len(hog_i) >= inferhog_min_hog_size_xml: # could be improved # hogs_a_rhog_xml = hog_i.to_orthoxml(**gene_id_name) hogs_a_rhog_xml_raw = hog_i.to_orthoxml(full_species_tree, score_flags) # @@ -135,6 +141,11 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder, pi if len(hog_i) > 1: logger.warning("issue 166312309 this is not a singleton "+str(hog_i._members)) + # any protein of the input rootHOG that is in no final HOG and was not reported as lost above + for rec in rhog_i: + if rec.id not in members_final: + ptrace.member_event(rec.id, "infer_subhogs", "missing_at_end", rhogid) + pickles_rhog_file = pickles_rhog_folder + '/file_' + rhogid + '.pickle' with open(pickles_rhog_file, 'wb') as handle: # dill_pickle.dump(hogs_rhogs_xml, handle, protocol=dill_pickle.HIGHEST_PROTOCOL) @@ -330,6 +341,7 @@ def __init__(self, node_species_tree:TreeNode, subhogs:List[HOG], rhogid:str, co self._base_msa_tree_filename = "HOG_"+rhogid+"_"+str(node_species_tree.name) self._outname_step = 0 self._removed_rep = set([]) + self._lost_members = set([]) # members dropped by split_hog at this level (for the protein trace) self._rep_lookup = self._prepare_lookups() self._msa_filter = self._instantiate_msa_filter() @@ -645,6 +657,12 @@ def merge_subhogs(self, reconciled_genetree:TreeNode, msa:MultipleSeqAlignment): split_parts = [list(n.name for n in sub.iter_leaves() if n.hogid == hogid) for sub in subtrees] split_hogs = split_hog(self.subhogs[hogid], self.node_species_tree.name, *split_parts) if split_hogs and len(split_hogs) > 1: + kept = set().union(*(h.get_members() for h in split_hogs)) + for member in self.subhogs[hogid].get_members() - kept: + self._lost_members.add(member) + ptrace.member_event(member, "infer_subhogs", "lost_in_split_hog", self.rhogid, + f"level={self.node_species_tree.name} hog={hogid} " + f"hog_size={len(self.subhogs[hogid])} parts={len(split_hogs)}") # we could split the current hog. self.subhogs.pop(hogid) self.subhogs.update({h.hogid: h for h in split_hogs}) @@ -673,6 +691,15 @@ def merge_subhogs(self, reconciled_genetree:TreeNode, msa:MultipleSeqAlignment): return new_hogs def process(self) -> List[HOG]: + members_before = set().union(set(), *(h.get_members() for h in list(self.subhogs.values()) + list(self._frozen_hogs_before))) + new_hogs = self._process() + members_after = set().union(set(), *(h.get_members() for h in new_hogs)) + for member in members_before - members_after - self._lost_members: + ptrace.member_event(member, "infer_subhogs", "lost_at_level", self.rhogid, + f"level={self.node_species_tree.name}") + return new_hogs + + def _process(self) -> List[HOG]: """Process the given roothog at the given taxonomic level. This method does all the processessing of the subhogs from the diff --git a/FastOMA/_protein_trace.py b/FastOMA/_protein_trace.py new file mode 100644 index 0000000..aae8c94 --- /dev/null +++ b/FastOMA/_protein_trace.py @@ -0,0 +1,82 @@ +"""Per-protein event log, so that the fate of every input protein can be traced. + +Each FastOMA step appends events to its own file in the task directory +(protein_events__.events, one per process so that parallel +Nextflow tasks never collide). The columns are, tab separated and without header: + + species protein step event roothog detail + +`fastoma-protein-trace` (FastOMA/protein_trace.py) merges all event files at the +end of the run into one row per input protein. + +Only changes are recorded (e.g. a protein being moved to another rootHOG, or +getting lost), plus the final state of every protein, so the files stay small. +""" + +import atexit +import threading +import uuid + +EVENT_FILE_SUFFIX = ".events" +COLUMNS = ["species", "protein", "step", "event", "roothog", "detail"] + +_lock = threading.Lock() +_handles = {} + + +def _handle(step): + fh = _handles.get(step) + if fh is None: + fh = open(f"protein_events_{step}_{uuid.uuid4().hex[:12]}{EVENT_FILE_SUFFIX}", "wt") + _handles[step] = fh + return fh + + +def event(species, protein, step, event_name, roothog="", detail=""): + line = "\t".join(str(x).replace("\t", " ") for x in (species, protein, step, event_name, roothog, detail)) + with _lock: + _handle(step).write(line + "\n") + + +def split_member(member): + """Internal FastOMA member ids look like 'protId||species||numericId', or several of + those joined by '_|_' for merged fragments. Returns a list of (species, protId).""" + out = [] + for part in str(member).split("_|_"): + fields = part.split("||") + if len(fields) >= 2: + out.append((fields[1], fields[0])) + return out + + +def member_event(member, step, event_name, roothog="", detail=""): + for species, protein in split_member(member): + event(species, protein, step, event_name, roothog, detail) + + +def rhog_assignment(rhogs_prots): + """(species, protId) -> rootHOG id, for diffing the rootHOG grouping steps.""" + return {sp_prot: rhogid for rhogid, sp_prots in rhogs_prots.items() for sp_prot in sp_prots} + + +def record_rhog_changes(before, after, step, event_prefix): + """Record proteins that changed rootHOG, or disappeared, between two assignments.""" + for sp_prot, old in before.items(): + new = after.get(sp_prot) + if new is None: + event(*sp_prot, step, f"{event_prefix}_removed", old) + elif new != old: + event(*sp_prot, step, f"{event_prefix}_moved", new, f"from={old}") + for sp_prot, new in after.items(): + if sp_prot not in before: + event(*sp_prot, step, f"{event_prefix}_added", new) + + +def close(): + with _lock: + for fh in _handles.values(): + fh.close() + _handles.clear() + + +atexit.register(close) diff --git a/FastOMA/_utils_roothog.py b/FastOMA/_utils_roothog.py index 222bc33..538c867 100644 --- a/FastOMA/_utils_roothog.py +++ b/FastOMA/_utils_roothog.py @@ -14,6 +14,7 @@ from .zoo.unionfind import UnionFind from . import logger +from . import _protein_trace as ptrace import collections filter_nonchild_rootHOG= False @@ -175,7 +176,12 @@ def parse_proteomes(folder=None, min_sequence_length=0): # list_oma_species #smallprot_recs_lists ={} for species_name in species_names: prot_address = os.path.join(folder, species_name + "." + fasta_format_keep) - prots_record = list(rec for rec in SeqIO.parse(prot_address, "fasta") if len(rec) >= min_sequence_length) + prots_record = [] + for rec in SeqIO.parse(prot_address, "fasta"): + if len(rec) >= min_sequence_length: + prots_record.append(rec) + else: + ptrace.event(species_name, rec.id, "infer_roothogs", "too_short", "", f"length={len(rec)} min={min_sequence_length}") #prots_record_small = list(rec for rec in SeqIO.parse(prot_address, "fasta") if len(rec) < min_sequence_length) # logger.debug(prots_record) logger.info(f"{species_name} contains {len(prots_record)} that are at least {min_sequence_length} long.") @@ -193,12 +199,23 @@ def create_root_hogs(hogmaps, conf): # Initial grouping rhogs_prots = group_prots_by_roothogs(hogmaps) logger.info(f"Initial grouping: {len(rhogs_prots)} root HOGs") + assignment = ptrace.rhog_assignment(rhogs_prots) + for (species_name, prot_id), rhogid in assignment.items(): + n_hits = len(hogmaps[species_name][prot_id]) + ptrace.event(species_name, prot_id, "infer_roothogs", "omamer_best_hit", rhogid, f"n_hits={n_hits}") # Handle singletons rhogs_prots = resolve_singletons(rhogs_prots, hogmaps, conf) + new_assignment = ptrace.rhog_assignment(rhogs_prots) + ptrace.record_rhog_changes(assignment, new_assignment, "infer_roothogs", "resolve_singletons") + assignment = new_assignment rhogs_prots = merge_rhogs2(hogmaps, rhogs_prots, conf) + new_assignment = ptrace.rhog_assignment(rhogs_prots) + ptrace.record_rhog_changes(assignment, new_assignment, "infer_roothogs", "merge_rhogs") + assignment = new_assignment rhogs_prots = filter_big_roothogs(hogmaps, rhogs_prots, conf) + ptrace.record_rhog_changes(assignment, ptrace.rhog_assignment(rhogs_prots), "infer_roothogs", "filter_big_roothogs") # # Merge similar HOGs # if conf.enable_merging: @@ -286,6 +303,7 @@ def parse_hogmap_omamer(proteomes, fasta_format_keep, folder=None): for row in reader: if row['hogid'] == "N/A": unmapped[species_name].append(row['qseqid']) + ptrace.event(species_name, row['qseqid'], "infer_roothogs", "no_omamer_hit") elif row['qseqid'] not in proteome: logger.debug(f"Ignoring {row} [{species_name}] from hogmap as not in proteome") else: @@ -511,8 +529,12 @@ def write_rhog(rhogs_prot_records, prot_recs_all, address_rhogs_folder, min_rhog # todo add the release id to file names rhogids_list[:2] > ['HOG:C0884658', 'HOG:C0709155'] SeqIO.write(rhog_recs, address_rhogs_folder + "/HOG_" + rhogid + ".fa", "fasta") rhogid_written_list.append(rhogid) + for rec in rhog_recs: + ptrace.member_event(rec.id, "infer_roothogs", "roothog_written", rhogid, f"size={len(rhog_recs)}") else: logger.debug("The roothog " +str(rhogid)+" was too small with size of "+str(len(rhog_recs))+" which is smaller than threshold "+str(min_rhog_size)) + for rec in rhog_recs: + ptrace.member_event(rec.id, "infer_roothogs", "roothog_too_small", rhogid, f"size={len(rhog_recs)} min={min_rhog_size}") #logger.debug(" issue 12314405 to check" ) #for prot1 in rhog_recs: #logger.debug("we are removing due to omamer signleton hog |*|" + str(prot1.id)) @@ -833,6 +855,8 @@ def collect_unmapped_singleton(rhogs_prots, unmapped, prot_recs_all, unmapped_si logger.debug(f"collected sequence records of {len(singleton_recs)} singleton proteins.") singleton_unmapped_recs = unmapped_recs + singleton_recs + for rec in singleton_unmapped_recs: + ptrace.member_event(rec.id, "infer_roothogs", "sent_to_clustering") SeqIO.write(singleton_unmapped_recs, unmapped_singleton_fasta , "fasta") logger.debug(f"Wrote {len(singleton_unmapped_recs)} sequences of unmapped and singleton proteins to {unmapped_singleton_fasta}.") return len(singleton_unmapped_recs) @@ -903,6 +927,7 @@ def write_clusters(fasta_clustered, address_rhogs_folder, min_rhog_size): for cluster in clusters: # cluster is like ['>sp|O67032|RF1_AQUAE', 'MLKEAYISRLDKLQEKYRKLQEELSKPEVIQD...', ...], i.e. an id and its sequence. # thus, check for double the min_rhog_size to have at least min_rhog_size proteins + member_ids = [cluster[pr_idi*2][1:].split(" ")[0] for pr_idi in range(int(len(cluster)/2))] if len(cluster) >= 2 * min_rhog_size: species_names = set([]) # it seems that genes from same species tend to cluster together, discard such clusters for pr_idi in range(int(len(cluster)/2)): @@ -912,7 +937,15 @@ def write_clusters(fasta_clustered, address_rhogs_folder, min_rhog_size): with open(address_rhogs_folder + "/HOG_clust" + str(cluster_iter) + ".fa", "w") as file_idx: for line in cluster: file_idx.write(line + "\n") + for member in member_ids: + ptrace.member_event(member, "infer_roothogs", "roothog_written", f"clust{cluster_iter}", f"size={len(member_ids)} source=linclust") cluster_iter += 1 + else: + for member in member_ids: + ptrace.member_event(member, "infer_roothogs", "cluster_single_species_discarded", "", f"size={len(member_ids)}") + else: + for member in member_ids: + ptrace.member_event(member, "infer_roothogs", "cluster_too_small", "", f"size={len(member_ids)}") return cluster_iter - 10000 diff --git a/FastOMA/protein_trace.py b/FastOMA/protein_trace.py new file mode 100644 index 0000000..0ae5549 --- /dev/null +++ b/FastOMA/protein_trace.py @@ -0,0 +1,216 @@ +"""Merge the per-step protein event files into one fate per input protein. + +fastoma-protein-trace --proteomes proteome --events events/ --roothog-tsv RootHOGs.tsv + +Outputs (in the current directory): + protein_fates.tsv one row per input protein, with exactly one final fate + protein_events.tsv.gz all recorded events, sorted by species and protein + protein_trace_summary.txt counts per fate, rootHOGs with most losses, consistency checks + +Fates, in the order they are decided: + too_short shorter than --min-sequence-length, never used + in_orthoxml member of a HOG of size >= 2 in the orthoXML + lost_in_split_hog dropped when split_hog did not cover all child subHOGs + lost_at_level dropped at a taxonomic level for another reason + final_singleton_not_written ended alone in its HOG; orthoXML only keeps HOGs of size >= 2 + missing_at_end in a rootHOG file, but in no final HOG and no loss recorded + roothog_not_processed in a rootHOG file, but no subHOG result (e.g. failed task) + singleton_roothog_dropped alone in its rootHOG and not placed in a cluster + no_omamer_hit_dropped no OMAmer hit and not placed in a cluster + removed_by_filter_big_roothogs removed from a very big rootHOG (low family score) + unaccounted no event at all (should not happen) +""" + +import argparse +import collections +import csv +import gzip +import os +import sys +from pathlib import Path + +from Bio import SeqIO + +from . import __version__ as fastoma_version +from ._protein_trace import COLUMNS, EVENT_FILE_SUFFIX + +LOSS_EVENTS = ("lost_in_split_hog", "lost_at_level") + + +def read_proteomes(folder): + """(species, protein) -> length, in input order.""" + proteins = {} + for f in sorted(os.listdir(folder)): + species, ext = f.rsplit(".", 1) + if ext not in ("fa", "fasta"): + continue + for rec in SeqIO.parse(os.path.join(folder, f), "fasta"): + proteins[(species, rec.id)] = len(rec) + return proteins + + +def read_events(paths): + events = collections.defaultdict(list) + n = 0 + for path in paths: + with open(path, "rt") as fh: + for line in fh: + fields = line.rstrip("\n").split("\t") + if len(fields) != len(COLUMNS): + continue + species, protein, step, event, roothog, detail = fields + events[(species, protein)].append((step, event, roothog, detail)) + n += 1 + return events, n + + +def event_files(inputs): + for p in inputs: + p = Path(p) + if p.is_dir(): + yield from sorted(p.rglob("*" + EVENT_FILE_SUFFIX)) + elif p.suffix == EVENT_FILE_SUFFIX: + yield p + + +def read_roothog_tsv(path): + """protein -> list of final RootHOG ids (a protein id could in principle occur in two species).""" + final = collections.defaultdict(list) + if path and os.path.exists(path): + with open(path, newline="") as fh: + for row in csv.DictReader(fh, delimiter="\t"): + final[row["Protein"]].append(row["RootHOG"]) + return final + + +def decide(evs): + """Return (fate, step, roothog, detail) from the events of one protein.""" + by_event = collections.defaultdict(list) + for step, event, roothog, detail in evs: + by_event[event].append((step, roothog, detail)) + + def first(name): + step, roothog, detail = by_event[name][0] + return name, step, roothog, detail + + if "too_short" in by_event: + return first("too_short") + if "in_final_hog" in by_event: + return first("in_final_hog") + for loss in LOSS_EVENTS: + if loss in by_event: + return first(loss) + if "final_singleton_not_written" in by_event: + return first("final_singleton_not_written") + if "missing_at_end" in by_event: + return first("missing_at_end") + if "roothog_written" in by_event: + _, step, roothog, detail = first("roothog_written") + return "roothog_not_processed", step, roothog, detail + if "filter_big_roothogs_removed" in by_event: + _, step, roothog, detail = first("filter_big_roothogs_removed") + return "removed_by_filter_big_roothogs", step, roothog, detail + reason = next((e for e in ("cluster_single_species_discarded", "cluster_too_small", "sent_to_clustering") if e in by_event), "") + if "no_omamer_hit" in by_event: + return "no_omamer_hit_dropped", "infer_roothogs", "", reason + if "roothog_too_small" in by_event or "omamer_best_hit" in by_event: + roothog = (by_event.get("roothog_too_small") or by_event["omamer_best_hit"])[0][1] + return "singleton_roothog_dropped", "infer_roothogs", roothog, reason + return "unaccounted", "", "", "" + + +def size_bucket(n): + if n is None: + return "-" + return "2-9" if n < 10 else "10-99" if n < 100 else "100-999" if n < 1000 else ">=1000" + + +def fastoma_protein_trace(): + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("--version", action="version", version="FastOMA v" + fastoma_version) + parser.add_argument("--proteomes", required=True, help="input proteome folder") + parser.add_argument("--events", required=True, nargs="+", help="event files or folders containing *.events") + parser.add_argument("--roothog-tsv", help="RootHOGs.tsv from fastoma-collect-subhogs (final RootHOG ids)") + parser.add_argument("--out-dir", default=".", help="where to write the three output files") + parser.add_argument("--strict", action="store_true", help="exit with an error if any protein is unaccounted") + conf = parser.parse_args() + + proteins = read_proteomes(conf.proteomes) + files = list(event_files(conf.events)) + events, n_events = read_events(files) + final_roothog = read_roothog_tsv(conf.roothog_tsv) + out = Path(conf.out_dir) + out.mkdir(parents=True, exist_ok=True) + + fates = collections.Counter() + fates_by_size = collections.Counter() + losses_by_rhog = collections.Counter() + rhog_size = {} + issues = collections.Counter() + + with open(out / "protein_fates.tsv", "w", newline="") as fh: + w = csv.writer(fh, delimiter="\t", lineterminator="\n") + w.writerow(["species", "protein", "length", "omamer_roothog", "omamer_n_hits", "roothog_changes", + "roothog", "roothog_size", "fate", "fate_step", "fate_detail", "final_roothog"]) + for key, length in proteins.items(): + evs = events.get(key, []) + fate, step, rhog, detail = decide(evs) + if fate == "in_final_hog": + fate = "in_orthoxml" + omamer = next(((r, d) for s, e, r, d in evs if e == "omamer_best_hit"), ("", "")) + changes = ";".join(f"{e.rsplit('_', 1)[0]}:{d.replace('from=', '')}->{r}" for s, e, r, d in evs if e.endswith("_moved")) + written = [(r, d) for s, e, r, d in evs if e == "roothog_written"] + if len(written) > 1: + issues["in_more_than_one_roothog_file"] += 1 + roothog, size = (written[0][0], int(written[0][1].split("size=")[1].split()[0])) if written else (rhog, None) + finals = final_roothog.get(key[1], []) if fate == "in_orthoxml" else [] + if fate == "in_orthoxml" and not finals: + issues["in_final_hog_but_not_in_RootHOGs.tsv"] += 1 + if len(finals) > 1: + issues["protein_id_in_several_RootHOGs.tsv_rows"] += 1 + w.writerow([key[0], key[1], length, omamer[0], omamer[1].replace("n_hits=", ""), changes, + roothog, "" if size is None else size, fate, step, detail, ",".join(finals)]) + fates[fate] += 1 + fates_by_size[(fate, size_bucket(size))] += 1 + if written: + rhog_size[roothog] = size + if fate in LOSS_EVENTS or fate in ("missing_at_end", "roothog_not_processed"): + losses_by_rhog[roothog] += 1 + + unknown = set(events) - set(proteins) + if unknown: + issues["events_for_proteins_not_in_proteomes"] = len(unknown) + + with gzip.open(out / "protein_events.tsv.gz", "wt") as fh: + fh.write("\t".join(COLUMNS) + "\n") + for key in sorted(events): + for ev in events[key]: + fh.write("\t".join((key[0], key[1]) + ev) + "\n") + + total = len(proteins) + buckets = ["2-9", "10-99", "100-999", ">=1000", "-"] + lines = [f"FastOMA protein trace ({len(files)} event files, {n_events:,} events)", "", + f"input proteins: {total:,}", "", "fate proteins %"] + for fate, n in fates.most_common(): + lines.append(f"{fate:32} {n:10,} {100 * n / total:8.2f}") + lines += ["", "fate by rootHOG size (proteins)", f"{'fate':32}" + "".join(f"{b:>10}" for b in buckets)] + for fate, _ in fates.most_common(): + lines.append(f"{fate:32}" + "".join(f"{fates_by_size[(fate, b)]:10,}" for b in buckets)) + lines += ["", "rootHOGs losing most proteins (lost, rootHOG size, rootHOG)"] + for rhog, n in losses_by_rhog.most_common(20): + lines.append(f" {n:7,} {rhog_size.get(rhog, 0):7,} {rhog}") + lines += ["", "consistency checks"] + lines.append(f" every input protein has exactly one fate: yes ({sum(fates.values()):,} rows)") + lines.append(f" unaccounted proteins: {fates.get('unaccounted', 0):,}") + for issue, n in sorted(issues.items()): + lines.append(f" {issue}: {n:,}") + summary = "\n".join(lines) + "\n" + (out / "protein_trace_summary.txt").write_text(summary) + print(summary) + + if conf.strict and fates.get("unaccounted", 0): + sys.exit(f"{fates['unaccounted']} proteins are unaccounted for") + + +if __name__ == "__main__": + fastoma_protein_trace() diff --git a/nextflow.config b/nextflow.config index 0efae03..ccec6b7 100644 --- a/nextflow.config +++ b/nextflow.config @@ -62,6 +62,7 @@ params { filter_gap_ratio_col = 0.5 nr_repr_per_hog = 5 min_sequence_length = 20 + protein_trace = false // write protein_fates.tsv: what happened to every input protein // scores written to FastOMA_HOGs.orthoxml -- set to disable a given score disable_completeness_score = false diff --git a/nextflow_schema.json b/nextflow_schema.json index 7be0e30..c77c16d 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -135,6 +135,11 @@ "min_sequence_length": { "type": "integer", "default": 20 + }, + "protein_trace": { + "type": "boolean", + "default": false, + "description": "Write protein_fates.tsv, protein_events.tsv.gz and protein_trace_summary.txt: what happened to every input protein." } } }, diff --git a/pyproject.toml b/pyproject.toml index ee77fb8..11af842 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,6 +44,7 @@ fastoma-collect-subhogs = "FastOMA.collect_subhogs:fastoma_collect_subhogs" fastoma-infer-roothogs = "FastOMA.infer_roothogs:fastoma_infer_roothogs" fastoma-infer-subhogs = "FastOMA.infer_subhogs:fastoma_infer_subhogs" fastoma-helper = "FastOMA.helper_scripts:main" +fastoma-protein-trace = "FastOMA.protein_trace:fastoma_protein_trace" [project.urls] Homepage = "https://github.com/DessimozLab/FastOMA" From d95c088524a583c8dd061cfa4234d7ac83e14158 Mon Sep 17 00:00:00 2001 From: aliyz Date: Mon, 19 Oct 2026 15:20:28 +0200 Subject: [PATCH 2/5] perf: protein trace only when enabled, stream the merge - Event logging in infer_roothogs and infer_subhogs is now off unless --protein-trace is given; Nextflow passes it when params.protein_trace is set. When off, no .events files are written and the extra bookkeeping (rootHOG assignment dicts, per-level member sets) is skipped. - fastoma-protein-trace sorts the event files on disk (GNU sort, stable) and processes one protein at a time instead of holding all events in memory; repeated strings (species, step, event, rootHOG) are interned. Peak memory on 500k synthetic proteins: 949 MB -> 224 MB, same output. protein_fates.tsv is now sorted by species and protein. - New withName "protein_trace" resources scaled by the event file size. --- FastOMA.nf | 9 ++- FastOMA/_infer_subhog.py | 31 +++++----- FastOMA/_protein_trace.py | 25 +++++++- FastOMA/infer_roothogs.py | 6 ++ FastOMA/infer_subhogs.py | 5 ++ FastOMA/protein_trace.py | 119 ++++++++++++++++++++++++++++---------- conf/base.config | 10 ++++ 7 files changed, 154 insertions(+), 51 deletions(-) diff --git a/FastOMA.nf b/FastOMA.nf index 5c3b6b3..2a44286 100644 --- a/FastOMA.nf +++ b/FastOMA.nf @@ -155,6 +155,7 @@ process infer_roothogs{ --splice ${splice_folder} \ --out-rhog-folder "omamer_rhogs" \ --min-sequence-length ${params.min_sequence_length} \ + ${ params.protein_trace.toString().toBoolean() ? "--protein-trace" : "" } \ -vv """ } @@ -227,7 +228,8 @@ process hog_big{ ${ params.write_genetrees ? "--gene-trees-write" : ""} \ ${ params.disable_completeness_score.toString().toBoolean() ? "--no-store-completeness-score" : "" } \ ${ params.disable_implied_losses_score.toString().toBoolean() ? "--no-store-implied-losses-score" : "" } \ - ${ params.disable_tcs_score.toString().toBoolean() ? "--no-store-tcs-score" : "" } + ${ params.disable_tcs_score.toString().toBoolean() ? "--no-store-tcs-score" : "" } \ + ${ params.protein_trace.toString().toBoolean() ? "--protein-trace" : "" } """ } @@ -264,7 +266,8 @@ process hog_rest{ ${ params.write_genetrees ? "--gene-trees-write" : ""} \ ${ params.disable_completeness_score.toString().toBoolean() ? "--no-store-completeness-score" : "" } \ ${ params.disable_implied_losses_score.toString().toBoolean() ? "--no-store-implied-losses-score" : "" } \ - ${ params.disable_tcs_score.toString().toBoolean() ? "--no-store-tcs-score" : "" } + ${ params.disable_tcs_score.toString().toBoolean() ? "--no-store-tcs-score" : "" } \ + ${ params.protein_trace.toString().toBoolean() ? "--protein-trace" : "" } """ } @@ -517,7 +520,7 @@ workflow { channel.empty().concat(pickle_big_rhog, pickle_rest_rhog).set{ all_rhog_pickle } (orthoxml_file, OrthologousGroupsFasta, OrthologousGroups_tsv, rootHOGs_tsv) = collect_subhogs(all_rhog_pickle.collect(), gene_id_dic_xml, omamer_rhogs, species_tree_checked, params.fasta_header_id_transformer) - if (params.protein_trace) { + if (params.protein_trace.toString().toBoolean()) { channel.empty() .concat(infer_roothogs.out.protein_events, hog_big.out.protein_events, hog_rest.out.protein_events) .flatten() diff --git a/FastOMA/_infer_subhog.py b/FastOMA/_infer_subhog.py index f49855e..8a2460b 100644 --- a/FastOMA/_infer_subhog.py +++ b/FastOMA/_infer_subhog.py @@ -108,10 +108,11 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder, pi members_final = set() for hog_i in hogs_a_rhog: tot_genes += len(hog_i) - members_final |= {part for m in hog_i.get_members() for part in m.split("_|_")} # merged fragments count per protein - final_event = "in_final_hog" if len(hog_i) >= inferhog_min_hog_size_xml else "final_singleton_not_written" - for member in hog_i.get_members(): - ptrace.member_event(member, "infer_subhogs", final_event, rhogid, f"hog={hog_i.hogid} size={len(hog_i)}") + if ptrace.is_enabled(): + members_final |= {part for m in hog_i.get_members() for part in m.split("_|_")} # merged fragments count per protein + final_event = "in_final_hog" if len(hog_i) >= inferhog_min_hog_size_xml else "final_singleton_not_written" + for member in hog_i.get_members(): + ptrace.member_event(member, "infer_subhogs", final_event, rhogid, f"hog={hog_i.hogid} size={len(hog_i)}") if len(hog_i) >= inferhog_min_hog_size_xml: # could be improved # hogs_a_rhog_xml = hog_i.to_orthoxml(**gene_id_name) hogs_a_rhog_xml_raw = hog_i.to_orthoxml(full_species_tree, score_flags) # @@ -142,9 +143,10 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder, pi logger.warning("issue 166312309 this is not a singleton "+str(hog_i._members)) # any protein of the input rootHOG that is in no final HOG and was not reported as lost above - for rec in rhog_i: - if rec.id not in members_final: - ptrace.member_event(rec.id, "infer_subhogs", "missing_at_end", rhogid) + if ptrace.is_enabled(): + for rec in rhog_i: + if rec.id not in members_final: + ptrace.member_event(rec.id, "infer_subhogs", "missing_at_end", rhogid) pickles_rhog_file = pickles_rhog_folder + '/file_' + rhogid + '.pickle' with open(pickles_rhog_file, 'wb') as handle: @@ -657,12 +659,13 @@ def merge_subhogs(self, reconciled_genetree:TreeNode, msa:MultipleSeqAlignment): split_parts = [list(n.name for n in sub.iter_leaves() if n.hogid == hogid) for sub in subtrees] split_hogs = split_hog(self.subhogs[hogid], self.node_species_tree.name, *split_parts) if split_hogs and len(split_hogs) > 1: - kept = set().union(*(h.get_members() for h in split_hogs)) - for member in self.subhogs[hogid].get_members() - kept: - self._lost_members.add(member) - ptrace.member_event(member, "infer_subhogs", "lost_in_split_hog", self.rhogid, - f"level={self.node_species_tree.name} hog={hogid} " - f"hog_size={len(self.subhogs[hogid])} parts={len(split_hogs)}") + if ptrace.is_enabled(): + kept = set().union(*(h.get_members() for h in split_hogs)) + for member in self.subhogs[hogid].get_members() - kept: + self._lost_members.add(member) + ptrace.member_event(member, "infer_subhogs", "lost_in_split_hog", self.rhogid, + f"level={self.node_species_tree.name} hog={hogid} " + f"hog_size={len(self.subhogs[hogid])} parts={len(split_hogs)}") # we could split the current hog. self.subhogs.pop(hogid) self.subhogs.update({h.hogid: h for h in split_hogs}) @@ -691,6 +694,8 @@ def merge_subhogs(self, reconciled_genetree:TreeNode, msa:MultipleSeqAlignment): return new_hogs def process(self) -> List[HOG]: + if not ptrace.is_enabled(): + return self._process() members_before = set().union(set(), *(h.get_members() for h in list(self.subhogs.values()) + list(self._frozen_hogs_before))) new_hogs = self._process() members_after = set().union(set(), *(h.get_members() for h in new_hogs)) diff --git a/FastOMA/_protein_trace.py b/FastOMA/_protein_trace.py index aae8c94..fe3521e 100644 --- a/FastOMA/_protein_trace.py +++ b/FastOMA/_protein_trace.py @@ -9,8 +9,9 @@ `fastoma-protein-trace` (FastOMA/protein_trace.py) merges all event files at the end of the run into one row per input protein. -Only changes are recorded (e.g. a protein being moved to another rootHOG, or -getting lost), plus the final state of every protein, so the files stay small. +Nothing is recorded unless enable() was called (the --protein-trace option of +fastoma-infer-roothogs and fastoma-infer-subhogs). Expect a few events per protein, +i.e. roughly 250-400 bytes per input protein in total. """ import atexit @@ -22,6 +23,16 @@ _lock = threading.Lock() _handles = {} +_enabled = False + + +def enable(): + global _enabled + _enabled = True + + +def is_enabled(): + return _enabled def _handle(step): @@ -33,6 +44,8 @@ def _handle(step): def event(species, protein, step, event_name, roothog="", detail=""): + if not _enabled: + return line = "\t".join(str(x).replace("\t", " ") for x in (species, protein, step, event_name, roothog, detail)) with _lock: _handle(step).write(line + "\n") @@ -50,17 +63,23 @@ def split_member(member): def member_event(member, step, event_name, roothog="", detail=""): + if not _enabled: + return for species, protein in split_member(member): event(species, protein, step, event_name, roothog, detail) def rhog_assignment(rhogs_prots): - """(species, protId) -> rootHOG id, for diffing the rootHOG grouping steps.""" + """(species, protId) -> rootHOG id, for diffing the rootHOG grouping steps (empty if tracing is off).""" + if not _enabled: + return {} return {sp_prot: rhogid for rhogid, sp_prots in rhogs_prots.items() for sp_prot in sp_prots} def record_rhog_changes(before, after, step, event_prefix): """Record proteins that changed rootHOG, or disappeared, between two assignments.""" + if not _enabled: + return for sp_prot, old in before.items(): new = after.get(sp_prot) if new is None: diff --git a/FastOMA/infer_roothogs.py b/FastOMA/infer_roothogs.py index 6d31fd2..a0fad31 100644 --- a/FastOMA/infer_roothogs.py +++ b/FastOMA/infer_roothogs.py @@ -3,6 +3,7 @@ from shutil import which from . import _utils_roothog, logger +from . import _protein_trace from . import __version__ as fastoma_version from .logging_setup import setup_logging @@ -34,9 +35,14 @@ def fastoma_infer_roothogs(): parser.add_argument("--big-rhog-size", required=False, type=int, default=50*1000, help= "For big rootHOGs, we have different heuristics") # big_rhog_size parser.add_argument("--big-fscore-thresh", required=False, type=int, default=95, help="For huge rootHOGs, we have different heuristics, like filtering low family score protiens") # big_fscore_thresh + parser.add_argument("--protein-trace", action="store_true", + help="write protein_events_*.events files, used by fastoma-protein-trace to trace every input protein") + conf = parser.parse_args() setup_logging(conf.v) logger.debug("Arguments: %s", conf) + if conf.protein_trace: + _protein_trace.enable() # Step 1: Parse input data logger.info("Parsing proteomes...") diff --git a/FastOMA/infer_subhogs.py b/FastOMA/infer_subhogs.py index 0445330..07608c9 100644 --- a/FastOMA/infer_subhogs.py +++ b/FastOMA/infer_subhogs.py @@ -3,6 +3,7 @@ from . import _utils_subhog from . import _infer_subhog +from . import _protein_trace from . import __version__ as fastoma_version from .logging_setup import setup_logging logger = logging.getLogger(__name__) @@ -55,11 +56,15 @@ def fastoma_infer_subhogs(): help="Store the ImpliedLosses score of each (Sub)HOG in the output orthoxml.") parser.add_argument("--store-tcs-score", action=argparse.BooleanOptionalAction, default=True, help="Store the TCSScore of each (Sub)HOG in the output orthoxml.") + parser.add_argument("--protein-trace", action="store_true", + help="write protein_events_*.events files, used by fastoma-protein-trace to trace every input protein") parser.add_argument('-v', action="count", default=0, help="Increase verbosity to info/debug") conf_infer_subhhogs = parser.parse_args() setup_logging(conf_infer_subhhogs.v, conf_infer_subhhogs.parallel) logger.debug("Arguments: %s", conf_infer_subhhogs) + if conf_infer_subhhogs.protein_trace: + _protein_trace.enable() address_rhogs_folder = conf_infer_subhhogs.input_rhog_folder # address_rhogs_folder = "./" # _config.input_rhog_folder diff --git a/FastOMA/protein_trace.py b/FastOMA/protein_trace.py index 0ae5549..3db536f 100644 --- a/FastOMA/protein_trace.py +++ b/FastOMA/protein_trace.py @@ -3,7 +3,7 @@ fastoma-protein-trace --proteomes proteome --events events/ --roothog-tsv RootHOGs.tsv Outputs (in the current directory): - protein_fates.tsv one row per input protein, with exactly one final fate + protein_fates.tsv one row per input protein, with exactly one final fate, sorted by species and protein protein_events.tsv.gz all recorded events, sorted by species and protein protein_trace_summary.txt counts per fate, rootHOGs with most losses, consistency checks @@ -19,14 +19,21 @@ no_omamer_hit_dropped no OMAmer hit and not placed in a cluster removed_by_filter_big_roothogs removed from a very big rootHOG (low family score) unaccounted no event at all (should not happen) + +The event files are sorted on disk (GNU sort), and the proteins are then processed +one at a time, so memory stays at a few hundred bytes per input protein no matter +how many events were recorded. """ import argparse import collections import csv import gzip +import itertools import os +import subprocess import sys +import tempfile from pathlib import Path from Bio import SeqIO @@ -38,30 +45,71 @@ def read_proteomes(folder): - """(species, protein) -> length, in input order.""" + """(species, protein) -> length.""" proteins = {} for f in sorted(os.listdir(folder)): species, ext = f.rsplit(".", 1) if ext not in ("fa", "fasta"): continue + species = sys.intern(species) for rec in SeqIO.parse(os.path.join(folder, f), "fasta"): proteins[(species, rec.id)] = len(rec) return proteins -def read_events(paths): - events = collections.defaultdict(list) - n = 0 - for path in paths: - with open(path, "rt") as fh: - for line in fh: - fields = line.rstrip("\n").split("\t") - if len(fields) != len(COLUMNS): - continue - species, protein, step, event, roothog, detail = fields - events[(species, protein)].append((step, event, roothog, detail)) - n += 1 - return events, n +def _parse_event_lines(lines): + _i = sys.intern # species, step, event and roothog repeat a lot, share one copy of each + for line in lines: + fields = line.rstrip("\n").split("\t") + if len(fields) != len(COLUMNS): + continue + species, protein, step, event, roothog, detail = fields + yield _i(species), protein, _i(step), _i(event), _i(roothog), detail + + +def sorted_events(paths, tmp_dir, sort_memory): + """Yield ((species, protein), [(step, event, roothog, detail), ...]), ordered by (species, protein). + + The event files are sorted on disk with GNU sort, so only the events of one protein + are held in memory at a time. The sort is stable and the files are given in name order + (infer_roothogs before infer_subhogs), so the events of a protein keep the order in + which they were written. LC_ALL=C makes sort compare bytes, which for UTF-8 is the + same order as comparing Python str tuples. + """ + if not paths: + return + with tempfile.NamedTemporaryFile("wb", dir=tmp_dir, suffix=".files0") as file_list: + file_list.write(b"".join(os.fsencode(p) + b"\0" for p in paths)) + file_list.flush() + cmd = ["sort", "-s", "-t", "\t", "-k1,1", "-k2,2", "-S", sort_memory, "-T", str(tmp_dir), + "--files0-from=" + file_list.name] + with subprocess.Popen(cmd, stdout=subprocess.PIPE, env=dict(os.environ, LC_ALL="C"), + encoding="utf-8") as proc: + for key, group in itertools.groupby(_parse_event_lines(proc.stdout), key=lambda f: (f[0], f[1])): + yield key, [f[2:] for f in group] + if proc.returncode: + raise RuntimeError(f"sorting the event files failed: {' '.join(cmd)} exited with {proc.returncode}") + + +def join_proteins_events(protein_keys, event_groups): + """Merge the sorted input proteins with the sorted event groups. + + Yields (key, events, is_input_protein); input proteins without any event get []. + """ + keys = iter(protein_keys) + key = next(keys, None) + for ev_key, evs in event_groups: + while key is not None and key < ev_key: + yield key, [], True + key = next(keys, None) + if key == ev_key: + yield key, evs, True + key = next(keys, None) + else: + yield ev_key, evs, False + while key is not None: + yield key, [], True + key = next(keys, None) def event_files(inputs): @@ -74,12 +122,12 @@ def event_files(inputs): def read_roothog_tsv(path): - """protein -> list of final RootHOG ids (a protein id could in principle occur in two species).""" - final = collections.defaultdict(list) + """protein -> tuple of final RootHOG ids (a protein id could in principle occur in two species).""" + final = {} if path and os.path.exists(path): with open(path, newline="") as fh: for row in csv.DictReader(fh, delimiter="\t"): - final[row["Protein"]].append(row["RootHOG"]) + final[row["Protein"]] = final.get(row["Protein"], ()) + (sys.intern(row["RootHOG"]),) return final @@ -133,14 +181,19 @@ def fastoma_protein_trace(): parser.add_argument("--roothog-tsv", help="RootHOGs.tsv from fastoma-collect-subhogs (final RootHOG ids)") parser.add_argument("--out-dir", default=".", help="where to write the three output files") parser.add_argument("--strict", action="store_true", help="exit with an error if any protein is unaccounted") + parser.add_argument("--tmp-dir", help="folder for the temporary files of the on-disk sort (default: --out-dir)") + parser.add_argument("--sort-memory", default="1G", help="main memory buffer of the on-disk sort (sort -S)") conf = parser.parse_args() proteins = read_proteomes(conf.proteomes) files = list(event_files(conf.events)) - events, n_events = read_events(files) final_roothog = read_roothog_tsv(conf.roothog_tsv) out = Path(conf.out_dir) out.mkdir(parents=True, exist_ok=True) + tmp_dir = Path(conf.tmp_dir) if conf.tmp_dir else out + tmp_dir.mkdir(parents=True, exist_ok=True) + n_events = 0 + n_unknown = 0 fates = collections.Counter() fates_by_size = collections.Counter() @@ -148,12 +201,21 @@ def fastoma_protein_trace(): rhog_size = {} issues = collections.Counter() - with open(out / "protein_fates.tsv", "w", newline="") as fh: + with open(out / "protein_fates.tsv", "w", newline="") as fh, \ + gzip.open(out / "protein_events.tsv.gz", "wt") as events_fh: w = csv.writer(fh, delimiter="\t", lineterminator="\n") w.writerow(["species", "protein", "length", "omamer_roothog", "omamer_n_hits", "roothog_changes", "roothog", "roothog_size", "fate", "fate_step", "fate_detail", "final_roothog"]) - for key, length in proteins.items(): - evs = events.get(key, []) + events_fh.write("\t".join(COLUMNS) + "\n") + event_groups = sorted_events(files, tmp_dir, conf.sort_memory) + for key, evs, is_input in join_proteins_events(sorted(proteins), event_groups): + for ev in evs: + events_fh.write("\t".join(key + ev) + "\n") + n_events += len(evs) + if not is_input: + n_unknown += 1 + continue + length = proteins[key] fate, step, rhog, detail = decide(evs) if fate == "in_final_hog": fate = "in_orthoxml" @@ -163,7 +225,7 @@ def fastoma_protein_trace(): if len(written) > 1: issues["in_more_than_one_roothog_file"] += 1 roothog, size = (written[0][0], int(written[0][1].split("size=")[1].split()[0])) if written else (rhog, None) - finals = final_roothog.get(key[1], []) if fate == "in_orthoxml" else [] + finals = final_roothog.get(key[1], ()) if fate == "in_orthoxml" else () if fate == "in_orthoxml" and not finals: issues["in_final_hog_but_not_in_RootHOGs.tsv"] += 1 if len(finals) > 1: @@ -177,15 +239,8 @@ def fastoma_protein_trace(): if fate in LOSS_EVENTS or fate in ("missing_at_end", "roothog_not_processed"): losses_by_rhog[roothog] += 1 - unknown = set(events) - set(proteins) - if unknown: - issues["events_for_proteins_not_in_proteomes"] = len(unknown) - - with gzip.open(out / "protein_events.tsv.gz", "wt") as fh: - fh.write("\t".join(COLUMNS) + "\n") - for key in sorted(events): - for ev in events[key]: - fh.write("\t".join((key[0], key[1]) + ev) + "\n") + if n_unknown: + issues["events_for_proteins_not_in_proteomes"] = n_unknown total = len(proteins) buckets = ["2-9", "10-99", "100-999", ">=1000", "-"] diff --git a/conf/base.config b/conf/base.config index f948658..69342a4 100644 --- a/conf/base.config +++ b/conf/base.config @@ -50,6 +50,16 @@ process { time = { 12.h * task.attempt * params.time_multiplier } } + withName: "protein_trace" { + cpus = { 1 } + memory = { + // ~400 B per protein are kept in memory, the event files have ~250-400 B per protein; + // the base also covers the 1 GB buffer of the on-disk sort + def events_bytes = [events].flatten().sum { f -> f.size() } ?: 0 + return (2.GB + events_bytes * 2) * task.attempt * params.memory_multiplier + } + } + withName: "collect_subhogs" { cpus = { 1 } } From fc9a45f56c272527c0502abeb5ff2e1076d0dd45 Mon Sep 17 00:00:00 2001 From: Ali Yazdizadeh Date: Wed, 21 Oct 2026 17:48:58 +0200 Subject: [PATCH 3/5] fix: add unit to the protein trace memory --- conf/base.config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/conf/base.config b/conf/base.config index 69342a4..6dd8fa9 100644 --- a/conf/base.config +++ b/conf/base.config @@ -56,7 +56,7 @@ process { // ~400 B per protein are kept in memory, the event files have ~250-400 B per protein; // the base also covers the 1 GB buffer of the on-disk sort def events_bytes = [events].flatten().sum { f -> f.size() } ?: 0 - return (2.GB + events_bytes * 2) * task.attempt * params.memory_multiplier + return (2.GB + (events_bytes * 2).B) * task.attempt * params.memory_multiplier } } From d27b4764743af984e07ad330e35acea513626f09 Mon Sep 17 00:00:00 2001 From: Ali Yazdizadeh Date: Wed, 21 Oct 2026 18:46:58 +0200 Subject: [PATCH 4/5] fix: log the not selected isoform in the traces --- FastOMA/_utils_roothog.py | 3 +++ FastOMA/protein_trace.py | 3 +++ 2 files changed, 6 insertions(+) diff --git a/FastOMA/_utils_roothog.py b/FastOMA/_utils_roothog.py index 538c867..f11dc65 100644 --- a/FastOMA/_utils_roothog.py +++ b/FastOMA/_utils_roothog.py @@ -1010,6 +1010,9 @@ def select_best_isoforms(species_names, isoform_by_gene_all, hogmaps): isoform_selected[species_name].append(protname_best) isoform_not_selected[species_name].extend([i for i in isoform_list if i != protname_best]) # flatten for all + for i in isoform_list: + if i != protname_best: + ptrace.event(species_name, i, "infer_roothogs", "isoform_not_selected", "", f"selected={protname_best}") return isoform_selected, isoform_not_selected diff --git a/FastOMA/protein_trace.py b/FastOMA/protein_trace.py index 3db536f..feb01a0 100644 --- a/FastOMA/protein_trace.py +++ b/FastOMA/protein_trace.py @@ -15,6 +15,7 @@ final_singleton_not_written ended alone in its HOG; orthoXML only keeps HOGs of size >= 2 missing_at_end in a rootHOG file, but in no final HOG and no loss recorded roothog_not_processed in a rootHOG file, but no subHOG result (e.g. failed task) + isoform_not_selected minor splicing variant, another isoform of the gene was selected singleton_roothog_dropped alone in its rootHOG and not placed in a cluster no_omamer_hit_dropped no OMAmer hit and not placed in a cluster removed_by_filter_big_roothogs removed from a very big rootHOG (low family score) @@ -155,6 +156,8 @@ def first(name): if "roothog_written" in by_event: _, step, roothog, detail = first("roothog_written") return "roothog_not_processed", step, roothog, detail + if "isoform_not_selected" in by_event: + return first("isoform_not_selected") if "filter_big_roothogs_removed" in by_event: _, step, roothog, detail = first("filter_big_roothogs_removed") return "removed_by_filter_big_roothogs", step, roothog, detail From a9dd74d6dae3e2d2a7390c6cc67d5d6e48c03088 Mon Sep 17 00:00:00 2001 From: Ali Yazdizadeh Date: Thu, 22 Oct 2026 11:05:32 +0200 Subject: [PATCH 5/5] fix: protein trace to use the same transformation as params.fasta_header_id_transformer --- FastOMA.nf | 6 ++++-- FastOMA/protein_trace.py | 7 ++++++- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/FastOMA.nf b/FastOMA.nf index 2a44286..2252849 100644 --- a/FastOMA.nf +++ b/FastOMA.nf @@ -313,6 +313,7 @@ process protein_trace { path proteome_folder path events, stageAs: "events/*" path roothogs_tsv + val id_transform output: path "protein_fates.tsv" @@ -323,7 +324,8 @@ process protein_trace { """ python -m FastOMA.protein_trace --proteomes ${proteome_folder} \ --events events/ \ - --roothog-tsv ${roothogs_tsv} + --roothog-tsv ${roothogs_tsv} \ + --id-transform ${id_transform} """ } @@ -526,7 +528,7 @@ workflow { .flatten() .collect() .set { all_protein_events } - protein_trace(proteome_folder, all_protein_events, rootHOGs_tsv) + protein_trace(proteome_folder, all_protein_events, rootHOGs_tsv, params.fasta_header_id_transformer) } c = hogmap.count().branch{ n-> TRUE: (n<=25 || params.force_pairwise_ortholog_generation) diff --git a/FastOMA/protein_trace.py b/FastOMA/protein_trace.py index feb01a0..0b79183 100644 --- a/FastOMA/protein_trace.py +++ b/FastOMA/protein_trace.py @@ -41,6 +41,7 @@ from . import __version__ as fastoma_version from ._protein_trace import COLUMNS, EVENT_FILE_SUFFIX +from .transformer import header_transformer LOSS_EVENTS = ("lost_in_split_hog", "lost_at_level") @@ -182,6 +183,8 @@ def fastoma_protein_trace(): parser.add_argument("--proteomes", required=True, help="input proteome folder") parser.add_argument("--events", required=True, nargs="+", help="event files or folders containing *.events") parser.add_argument("--roothog-tsv", help="RootHOGs.tsv from fastoma-collect-subhogs (final RootHOG ids)") + parser.add_argument("--id-transform", choices=("UniProt", "noop"), default="noop", + help="ID transformer that fastoma-collect-subhogs used for the protein IDs in RootHOGs.tsv") parser.add_argument("--out-dir", default=".", help="where to write the three output files") parser.add_argument("--strict", action="store_true", help="exit with an error if any protein is unaccounted") parser.add_argument("--tmp-dir", help="folder for the temporary files of the on-disk sort (default: --out-dir)") @@ -191,6 +194,7 @@ def fastoma_protein_trace(): proteins = read_proteomes(conf.proteomes) files = list(event_files(conf.events)) final_roothog = read_roothog_tsv(conf.roothog_tsv) + id_transformer = header_transformer(conf.id_transform) out = Path(conf.out_dir) out.mkdir(parents=True, exist_ok=True) tmp_dir = Path(conf.tmp_dir) if conf.tmp_dir else out @@ -228,7 +232,8 @@ def fastoma_protein_trace(): if len(written) > 1: issues["in_more_than_one_roothog_file"] += 1 roothog, size = (written[0][0], int(written[0][1].split("size=")[1].split()[0])) if written else (rhog, None) - finals = final_roothog.get(key[1], ()) if fate == "in_orthoxml" else () + # RootHOGs.tsv has the protein IDs as written to the orthoXML, i.e. after the id transform + finals = final_roothog.get(id_transformer.transform(key[1]), ()) if fate == "in_orthoxml" else () if fate == "in_orthoxml" and not finals: issues["in_final_hog_but_not_in_RootHOGs.tsv"] += 1 if len(finals) > 1: