Tcs and implied loss scores

FastOMA.nf

@@ -222,7 +222,10 @@ process hog_big{
                                --gap-ratio-col ${params.filter_gap_ratio_col} \
                                --number-of-samples-per-hog ${params.nr_repr_per_hog} \
                                ${ params.write_msas ? "--msa-write" : ""} \
-                               ${ params.write_genetrees ? "--gene-trees-write" : ""}
+                               ${ 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" : "" }
     """
 }
 
@@ -255,7 +258,10 @@ process hog_rest{
                               --gap-ratio-col ${params.filter_gap_ratio_col} \
                               --number-of-samples-per-hog ${params.nr_repr_per_hog} \
                               ${ params.write_msas ? "--msa-write" : ""} \
-                              ${ params.write_genetrees ? "--gene-trees-write" : ""}
+                              ${ 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" : "" }
     """
 }
 

FastOMA/_infer_subhog.py

@@ -26,7 +26,7 @@
 from . import _wrappers, logger
 from . import _utils_subhog
 from . import _utils_frag_SO_detection
-from ._hog_class import HOG, Representative, split_hog
+from ._hog_class import HOG, Representative, split_hog, attach_scores, _member_species, _species_name_index, ScoreFlags
 from ._utils_subhog import MSAFilter, MSAFilterElbow, MSAFilterTrimAL
 
 from .zoo.utils import unique
@@ -73,7 +73,7 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder,  pi
     # the file "species_tree_checked.nwk" is created by the check_input.py
     (species_tree) = _utils_subhog.read_species_tree(conf_infer_subhhogs.species_tree)
 
-    (species_tree, species_names_rhog, prot_names_rhog) = _utils_subhog.prepare_species_tree(rhog_i, species_tree, rhogid)
+    (species_tree, species_names_rhog, prot_names_rhog, full_species_tree) = _utils_subhog.prepare_species_tree(rhog_i, species_tree, rhogid)
     species_names_rhog = list(set(species_names_rhog))
     logger.info("Number of unique species in rHOG " + rhogid + " is " + str(len(species_names_rhog)) + ".")
 
@@ -96,20 +96,27 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder,  pi
     if not keep_subhog_each_pickle:
         shutil.rmtree(pickles_subhog_folder)
 
+    score_flags = ScoreFlags(
+        store_completeness_score=conf_infer_subhhogs.store_completeness_score,
+        store_implied_losses_score=conf_infer_subhhogs.store_implied_losses_score,
+        store_tcs_score=conf_infer_subhhogs.store_tcs_score,
+    )
+
     tot_genes, placed_genes = 0, 0
     hogs_rhogs_xml = []
     for hog_i in hogs_a_rhog:
         tot_genes += 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()    # <generef  >      <paralg object >
+            hogs_a_rhog_xml_raw = hog_i.to_orthoxml(full_species_tree, score_flags)    # <generef  >      <paralg object >
             if orthoxml_v03 and 'paralogGroup' in str(hogs_a_rhog_xml_raw) :
                 # in version v0.3 of orthoxml, there shouldn't be any paralogGroup at root level. Let's put them inside an orthogroup should be in
                 hog_elemnt = ET.Element('orthologGroup', attrib={"id": str(hog_i.hogid)})
-                num_species_tax_hog = len(set([i.split("||")[1] for i in hog_i.get_members()]))
-                completeness_score = round(num_species_tax_hog / hog_i.taxlevel.size, 4)
-                ET.SubElement(hog_elemnt, "score",
-                              attrib={"id": "CompletenessScore", "value": str(completeness_score)})
+                species_of_members = _member_species(hog_i.get_members())
+                scoring_node = hog_i.taxlevel
+                if full_species_tree is not None:
+                    scoring_node = _species_name_index(full_species_tree)[hog_i.taxlevel.name]
+                attach_scores(hog_elemnt, hog_i, scoring_node, species_of_members, score_flags=score_flags)
                 ET.SubElement(hog_elemnt, "property", attrib={"name": "TaxRange", "value": str(hog_i.taxname)})
                 hog_elemnt.append(hogs_a_rhog_xml_raw)
                 hogs_a_rhog_xml = hog_elemnt
@@ -163,6 +170,10 @@ def build_xml_from_rhog(rhogid:str, seqs:List[SeqRecord], hogs:List[ET.Element])
     scores = ET.SubElement(root, "scores")
     ET.SubElement(scores, "scoreDef",
                   {"id": "CompletenessScore","desc": "Fraction of expected species with genes in the (Sub)HOG"})
+    ET.SubElement(scores, "scoreDef",
+                  {"id": "TCSScore", "desc": "Taxonomic Congruence Score: how well the (Sub)HOG structure matches the species tree topology"})
+    ET.SubElement(scores, "scoreDef",
+                  {"id": "ImpliedLosses", "desc": "Number of implied gene loss events (Dollo parsimony) within the (Sub)HOG's taxonomic range"})
     groups = ET.SubElement(root, 'groups')
     for hog in hogs:
         groups.append(hog)

FastOMA/_utils_subhog.py

@@ -253,7 +253,10 @@ def prepare_species_tree(rhog_i: List[SeqRecord], species_tree: Tree, rhogid: st
     orthoxml_to_newick.py function for extracting orthoxml_to_newick.py subtree from the input species tree  orthoxml_to_newick.py.k.orthoxml_to_newick.py pruning,
     based on the names of species in the rootHOG.
 
-    output: species_tree (pruned), species_names_rhog, prot_names_rhog
+    output: species_tree (pruned), species_names_rhog, prot_names_rhog, full_species_tree (unpruned,
+    with the `size` feature set) -- kept around so that scores computed later (e.g. ImpliedLosses,
+    TCSScore) can be evaluated against the real species tree topology, not just the species present
+    in this rootHOG.
     """
     assert len(rhog_i) > 0, 'input hog_i is empty, probably previous step find_rhog has issue, rhogs/HOG_B0'+rhogid+'is empty?'
     species_names_rhog = []
@@ -275,9 +278,11 @@ def prepare_species_tree(rhog_i: List[SeqRecord], species_tree: Tree, rhogid: st
     for n in species_tree.traverse():
         n.add_feature("size", len(n))
 
+    full_species_tree = species_tree.copy()
+
     mrca = species_tree.get_common_ancestor(species_names_uniqe)
     mrca.prune(species_names_uniqe, preserve_branch_length=True)
-    return mrca, species_names_rhog, prot_names_rhog
+    return mrca, species_names_rhog, prot_names_rhog, full_species_tree
 
 
 def label_sd_internal_nodes(tree_out, threshold_dubious_sd):

FastOMA/collect_subhogs.py

@@ -26,6 +26,24 @@
 
 # This code collect subhogs and writes outputs.
 
+SCORE_DEFS = {
+    "CompletenessScore": "Fraction of expected species with genes in the (Sub)HOG",
+    "TCSScore": "Taxonomic Congruence Score: how well the (Sub)HOG structure matches the species tree topology",
+    "ImpliedLosses": "Number of implied gene loss events (Dollo parsimony) within the (Sub)HOG's taxonomic range",
+}
+
+
+def peek_score_ids(pickle_folder: Path) -> set:
+    """Returns the set of score ids attached to the first HOG found in pickle_folder.
+
+    Score ids depend only on the (uniform, run-wide) --store-*-score flags passed to
+    fastoma-infer-subhogs, not on individual HOGs, so a single HOG is representative of
+    the whole run and stops the scan after the first non-empty pickle file."""
+    for hog in iter_hogs(pickle_folder):
+        return {child.get('id') for child in hog if child.tag == 'score'}
+    return set()
+
+
 def iter_hogs(pickle_folder: Path):
     cnt = 0
     nr_hogs = 0
@@ -102,7 +120,14 @@ def _annotateGroupR(node: ET.ElementTree, og: str, idx: int = 0):
 
     omamer_roothog_id = ":".join(hog.get('id').split("_")[0:2])
     fam_elem = ET.Element("property", {"name": "OMAmerRootHOG", "value": omamer_roothog_id})
-    hog.insert(1, fam_elem)
+    # orthoxml requires all <score> children before any <property> children, so insert
+    # the new property right after the trailing run of <score> elements, not at a fixed index.
+    insert_idx = 0
+    for child in hog:
+        if child.tag != "score":
+            break
+        insert_idx += 1
+    hog.insert(insert_idx, fam_elem)
     _annotateGroupR(hog, "HOG:{:07d}".format(fam))
     return hog
 
@@ -186,9 +211,12 @@ def write_hog_orthoxml(pickle_folder, output_xml_name, gene_id_pickle_file, id_t
     logger.debug("gene_xml is created.")
     orthoxml_file.append(taxonomy)
 
-    scores = ET.SubElement(orthoxml_file, "scores")
-    ET.SubElement(scores, "scoreDef", {"id": "CompletenessScore",
-                                       "desc": "Fraction of expected species with genes in the (Sub)HOG"})
+    present_score_ids = peek_score_ids(Path(pickle_folder))
+    if present_score_ids:
+        scores = ET.SubElement(orthoxml_file, "scores")
+        for score_id, desc in SCORE_DEFS.items():
+            if score_id in present_score_ids:
+                ET.SubElement(scores, "scoreDef", {"id": score_id, "desc": desc})
 
     #  #### create the groups of orthoxml   ####
     groups_xml = ET.SubElement(orthoxml_file, "groups")

FastOMA/infer_subhogs.py

@@ -49,6 +49,12 @@ def fastoma_infer_subhogs():
                         help="For trimming the MSA, the threshold of ratio of gaps for each column.")
     parser.add_argument("--min-col-trim", required=False, type=int, default=50,  # todo min rows trim
                         help="min no. columns in msa to consider for filtering")
+    parser.add_argument("--store-completeness-score", action=argparse.BooleanOptionalAction, default=True,
+                        help="Store the CompletenessScore of each (Sub)HOG in the output orthoxml.")
+    parser.add_argument("--store-implied-losses-score", action=argparse.BooleanOptionalAction, default=True,
+                        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('-v', action="count", default=0, help="Increase verbosity to info/debug")
     conf_infer_subhhogs = parser.parse_args()
 

FastOMA/zoo/wrappers/treebuilders/fasttree.py

@@ -94,7 +94,7 @@ def _call(self, filename, *args, **kwargs):
             self.stderr = self.cli.get_stderr()
             last_error_line = self.stderr.split('\n')[-1].strip()
             logger.error('FastTree returned non-zero exit status: {}'.format(self.returncode))
-            logger.error('Output of FastTree:\n\n%s\nstdout=\n%s\n{}\n\n%s\nstderr=\n%s\n{}\n\n',
+            logger.error('Output of FastTree:\n\n%s\nstdout=\n%s\n%s\n\n%s\nstderr=\n%s\n%s\n\n',
                          "=" * 30, "=" * 30, summarize_long_message(self.stdout),
                          "=" * 30, "=" * 30, summarize_long_message(self.stderr))
             if self.returncode < 0:

FastOMA/zoo/wrappers/treebuilders/parsers.py

@@ -7,9 +7,8 @@
 logger.addHandler(logging.StreamHandler())
 
 
-
-FLOAT = Word(nums + '.-').setParseAction(lambda x: float(x[0]))
-INT = Word(nums).setParseAction(lambda x: int(x[0]))
+FLOAT = Word(nums + '.-').set_parse_action(lambda x: float(x[0]))
+INT = Word(nums).set_parse_action(lambda x: int(x[0]))
 WORD = Word(alphanums + '_')
 SPACEDWORD = Word(alphanums + ' _')
 MODEL_CONCAT_WORD = Word(alphanums + "+_-")

conf/test-fungi.config

@@ -1,7 +1,7 @@
 // Default configuration for Nextflow
 
 params {
-    test_data_url         = "https://zenodo.org/records/17434495/files/fungi-30.tgz?download=1"
+    input                 = "https://zenodo.org/records/17434495/files/fungi-30.tgz?download=1"
     report                = true
-    omamer_db             = "${projectDir}/testdata/test.h5"
+    omamer_db             = "https://zenodo.org/records/20814376/files/saccharomyceta.h5?download=1"
 }

conf/test-mammalia.config

@@ -1,6 +1,7 @@
 // Default configuration for Nextflow
 
 params {
-    test_data_url         = "https://zenodo.org/records/17434495/files/mammalia-22.tgz?download=1"
+    input                 = "https://zenodo.org/records/17434495/files/mammalia-22.tgz?download=1"
+    omamer_db             = "https://zenodo.org/records/20814376/files/Metazoa.h5?download=1"
     report                = true
 }

nextflow.config

@@ -63,6 +63,10 @@ params {
   nr_repr_per_hog = 5
   min_sequence_length = 20
 
+  // scores written to FastOMA_HOGs.orthoxml -- set to disable a given score
+  disable_completeness_score    = false
+  disable_implied_losses_score  = false
+  disable_tcs_score             = false
 
   // other parameters
   debug_enabled = false   // generates additional files for debugging

nextflow_schema.json

@@ -171,6 +171,18 @@
           "type": "boolean",
           "description": "Force generation of pairwise orthologs even for large datasets (may be slow)",
           "help": "By default, pairwise orthologs are not generated for datasets with more than 25 species to avoid long runtimes. Enable this option to override this behavior."
+        },
+        "disable_completeness_score": {
+          "type": "boolean",
+          "description": "Disable computing/storing the CompletenessScore of each (Sub)HOG in the output orthoxml"
+        },
+        "disable_implied_losses_score": {
+          "type": "boolean",
+          "description": "Disable computing/storing the ImpliedLosses score of each (Sub)HOG in the output orthoxml"
+        },
+        "disable_tcs_score": {
+          "type": "boolean",
+          "description": "Disable computing/storing the TCSScore of each (Sub)HOG in the output orthoxml"
         }
       }
     },

nf-tests/default.nf.test.snap

@@ -50,7 +50,7 @@
             [
                 {
                     "name": "FastOMA_HOGs.orthoxml",
-                    "lineCount": 162
+                    "lineCount": 194
                 },
                 {
                     "name": "OrthologousGroups.tsv",
@@ -80,8 +80,8 @@
         ],
         "meta": {
             "nf-test": "0.9.3",
-            "nextflow": "24.10.5"
+            "nextflow": "26.04.6"
         },
-        "timestamp": "2025-10-23T12:09:08.301966"
+        "timestamp": "2026-08-25T18:07:32.287389"
     }
 }

tests/data/tcs_taxonomy.tsv

@@ -0,0 +1,5 @@
+query	lineage
+SP1	1 (class), 2 (order), 3 (family), 4 (genus), 5 (species)
+SP2	1 (class), 2 (order), 3 (family), 4 (genus), 6 (species)
+SP3	1 (class), 7 (order), 8 (family), 9 (genus), 10 (species)
+SP4	1 (class), 7 (order), 11 (family), 12 (subfamily), 13 (genus), 14 (species), 15 (subspecies)

tests/data/tcs_tree1.nwk

@@ -0,0 +1 @@
+((SP1,SP2),(SP3,SP4));

tests/data/tcs_tree2.nwk

@@ -0,0 +1 @@
+(((SP1,SP2),SP3),SP4);

tests/data/tcs_tree3.nwk

@@ -0,0 +1 @@
+((SP1,SP3),(SP2,SP4));