bio-comparative-genomics-ortholog-inference

## Version Compatibility Reference examples tested with: OrthoFinder 3.0+ (Emms et al 2026 Nat Methods 23:1327), SonicParanoid 2.0.8+ (Cosentino 2024), Broccoli 1.2+ (Derelle 2020), ProteinOrtho 6.3.0+ (Lechner 2011 + recent), OMA standalone 2.6.0+, FastOMA 0.3.5+ (Majidian 2025), eggNOG-mapper 2.1.12+, JustOrthologs 2.0+, DIAMOND 2.1.10+, MMseqs2 17-b804f+, IQ-TREE 2.3.6+, BUSCO 5.7+, Compleasm 0.2.7+, BioPython 1.84+, R 4.4+ for downstream tree-based reconciliation. Before using code patterns, verify installed versions match. If versions differ: - CLI: `orthofinder --help`; `sonicparanoid --help`; `oma --help` - Python: `pip show eggnog-mapper`; `which fastoma` If code throws `Diamond requires N more sequences than provided`, `KeyError on species tree taxa`, `STAG branch length 0`, or `HOG file format mismatch`, the OrthoFinder v2 -> v3 file layout changed (Orthogroups/ -> Phylogenetic_Hierarchical_Orthogroups/; rooted gene trees are now per-HOG); update parsing accordingly. # Ortholog Inference **"Find the orthologs of my gene(s) across these species"** -> Choose between graph-based (RBH / similarity-clustering: fast, lower recall) and tree-based (gene-tree reconciliation: higher accuracy, slower) frameworks; recognize that "orthology" splits into 1:1, 1:many, many:many, and the practical unit for most pipelines is the **HOG (Hierarchical Orthologous Group)** -- a maximal cluster of genes descended from a single ancestral gene at a defined taxonomic level (Altenhoff 2013 PLoS ONE 8:e53786). The "ortholog conjecture" (orthologs more functionally similar than paralogs) is supported but weakly (Altenhoff 2012 PLoS Comp Biol 8:e1002514); don't treat 1:1 ortholog labeling as automatic functional equivalence. - CLI: `orthofinder -f proteomes/ -t 16 -M msa` -- HOG output in v3 layout - CLI: `sonicparanoid -i proteomes/ -o output --mode default` -- ML predictor + protein language model - CLI: `broccoli.py -dir proteomes/ -threads 16` -- direct OG with chimeric handling - CLI: `oma standalone` HOG inference at every taxonomic level - CLI: `proteinortho6.pl --project=run proteomes/*.faa` -- graph clustering with optional synteny - CLI: `emapper.py -i proteins.faa --output project --cpu 16` -- eggNOG annotation transfer ## Algorithmic Taxonomy | Tool | Approach | Output | Strength | Fails when | |------|----------|--------|----------|------------| | OrthoFinder3 (Emms et al 2026 Nat Methods) | DIAMOND2-ultra search -> gene trees -> rooted via STRIDE -> HOG inference at every node | HOGs at every taxonomic level, orthologs, species tree, gene-duplication events | Best Quest-for-Orthologs benchmark (Altenhoff 2024); HOGs are ~7% more accurate than v2 orthogroups (5-7% higher f-score); outgroup inclusion further improves accuracy | Slow with > 200 species without `-c 1` clustering pre-step; new file layout breaks old parsers | | SonicParanoid2 (Cosentino 2024 GB 25:195) | RBH + ML classifier with protein language model embeddings | Pairwise orthologs + orthogroups | Best accuracy among graph methods at competitive speed; ML correction reduces InParanoid errors | InParanoid lineage still has paralogy confusion in WGD clades; close-relative duplications hard to resolve | | Broccoli (Derelle 2020 MBE 37:3389) | k-mer similarity -> directed graph -> chimera-aware OG inference | OGs + chimera flagging | Robust to chimeric assemblies; runs without species tree | Less accurate than OrthoFinder on benchmark; no HOG output | | ProteinOrtho 6 (Lechner 2011 BMC Bioinf 12:124) | Pairwise BLAST/DIAMOND + connectivity graph; optional synteny module | Orthogroups + synteny option | Fast; scales to 1000+ genomes; `-synteny` enables co-linear-anchor filtering | Lower recall than tree-based; synteny module slow and requires GFFs | | OMA standalone (Altenhoff 2019 GR 29:1152) | Strict RBH + verification + HOG inference | HOG database; orthologs at each taxonomic level | Conservative; highest precision in QfO benchmarks; "Fast" mode for prefiltering | Lowest recall among methods (Altenhoff 2024); slow for large datasets | | FastOMA (Majidian 2025 Nat Methods 22:269) | OMA HOG inference with GPU-accelerated DIAMOND + Roothap | Same HOG output as OMA, 10-100x faster | Scales OMA to 1000+ genomes | Newer; less benchmarked in production | | eggNOG-mapper 2 (Cantalapiedra 2021 MBE 38:5825) | DIAMOND/MMseqs2 against eggNOG 5 -> map to precomputed orthogroups | OGs + functional annotation (GO/KEGG/COG) | Standard for functional annotation propagation; phylogeny-aware | Pre-computed OGs; cannot add novel species coherently; only as fresh as eggNOG release | | JustOrthologs (Miller 2019 Bioinformatics 35:546) | DNA-based; exon-aware; close-species RBH | Pairwise orthologs | Extremely fast for closely related species (same family); preserves splice variants | Only suitable for closely related species | | TOGA (Kirilenko 2023 Science 380:eabn3107) | Whole-genome-alignment chain -> ML projection + intactness classification | Per-query orthologs with intact/lost/missing call | Modern paradigm for vertebrate-scale orthology; handles gene loss explicitly; integrates with CESAR 2.0 | Requires WGA (Cactus); not designed for prokaryotes or fungi | | HOGENOM / HOGsuite (Penel 2009 BMC Bioinf 10:S3) | Tree-based HOGs in databases | Pre-computed HOG database | Legacy; for downstream use of stored HOGs | Not for new computation; outdated taxon sampling | Methodology evolves; Quest-for-Orthologs benchmarks (Altenhoff 2024 NAR Genom Bioinform 6:lqae167) refresh annually. Verify the current QfO benchmark results before locking on a single tool for novel benchmarking-grade work. ## Decision Tree by Experimental Scenario | Scenario | Recommended approach | Why | |----------|------------------------|-----| | 5-50 vertebrate / animal genomes for phylogenomics | OrthoFinder3 `-M msa -A mafft -T iqtree` | HOG quality at moderate scale; species tree included | | 50-200 genomes, any clade | OrthoFinder3 default; SonicParanoid2 as second method | Cross-validate consensus HOGs | | 200-1000 genomes, scaling required | FastOMA or SonicParanoid2 | OrthoFinder3 slow; cluster-first option (`-c 1`) helps | | 1000+ bacterial genomes (pangenome scope) | Use [[pangenome-analysis]] pipeline; OrthoFinder unsuitable | Pangenome-specific tools needed | | Closely related strains / same family | JustOrthologs or ProteinOrtho with synteny | Fast; preserves splice variants | | WGD-affected lineage (plants, salmonids, teleosts) | OrthoFinder3 + synteny verification (see [[synteny-analysis]]); or GENESPACE | WGD inflates paralogy; synteny anchoring required | | Functional annotation transfer | eggNOG-mapper 2 | Pre-computed; phylogeny-aware; integrates GO/KEGG | | Single-copy orthologs for concatenation phylogenomics | OrthoFinder3 `Phylogenetic_Hierarchical_Orthogroups/N0.tsv` filter for 1:1 | Most rigorous HOG layer; standard concatenation input | | Gene-loss detection across mammals/birds at scale | TOGA + CESAR 2.0 | Explicit intact/lost classification; handles assembly-gap noise | | Functional annotation transfer in poorly characterized genome | eggNOG-mapper + OrthoFinder3 single-copy orthologs cross-validated | Functional priors + phylogenetic confidence | | Ortholog detection in highly fragmented assembly | BUSCO/Compleasm first to assess completeness; flag affected OGs | Assembly fragmentation creates false absence -> spurious lineage-specific losses | | Distant homolog detection (50-200 Myr divergence) | MMseqs2 sensitive (`-s 7.5`) + OrthoFinder3 | Default DIAMOND misses distant homologs; sensitive search needed | | Phylogenetic orthology with paralog-tolerant species trees | ASTRAL-Pro2 on OrthoFinder gene trees | Coalescent species tree from gene trees; handles paralogy explicitly | | HOG-level functional propagation across many species | OMA standalone HOG database | Strict hierarchical orthology; supports functional propagation per taxonomic level | | Synteny-anchored orthology in repeat-heavy genome | ProteinOrtho `-synteny` or GENESPACE | Filters tandem duplicates; uses gene order to disambiguate paralogs | ## Per-Method Failure Modes ### Hidden paralogy from missing outgroup **Trigger:** OrthoFinder run without a sufficiently distant outgroup; only ingroup taxa included. **Mechanism:** OrthoFinder roots gene trees via STRIDE / STAG using outgroup-based duplication signals. Without an outgroup, the root inferred for each gene tree may be internal, mistakenly classifying a duplication-then-loss pattern as orthology. The "1:1 orthologs" returned may actually be hidden paralogs (Emms & Kelly 2017 MBE 34:3267 STRIDE). **Symptom:** Phylogenomic concatenation produces low-bootstrap species tree; per-gene trees show inconsistent rooting; same-clade species show longer-than-expected branches in single-copy ortholog trees. **Fix:** Include >= 1 outgroup taxon at the next-higher taxonomic level (sister phylum / class / order). For vertebrates, use cyclostomes (lamprey) as outgroup for jawed vertebrates; for plants, use a non-flowering plant for angiosperms. Re-run with outgroup; cross-check HOG file `Phylogenetic_Hierarchical_Orthogroups/N0.tsv` for 1:1 stability. ### Splice isoforms inflating copy number **Trigger:** Proteome FASTA contains multiple isoforms per gene (Ensembl, RefSeq with `-NR` flag). **Mechanism:** OrthoFinder / SonicParanoid / OMA treat each protein sequence as a gene; isoforms are clustered together within orthogroups but inflate the apparent gene count per species, producing artifactual "co-orthologs" that are just isoforms of one gene. **Symptom:** Per-species gene count > 2x what gene-annotation pipeline reported; orthogroups contain multiple proteins from same gene; gene names contain isoform suffixes (`.1`, `-iso1`, etc.). **Fix:** Pre-filter proteomes to longest isoform per gene. Use `agat_sp_keep_longest_isoform.pl` (AGAT toolkit), Biopython snippet on Ensembl gene-isoform mapping, or `Trinotate` longest-ORF picker. OrthoFinder3 ships a wrapper `tools/primary_transcript.py` for Ensembl/UniProt format. ### Annotation heterogeneity inflating lineage-specific OGs **Trigger:** Genomes annotated by different pipelines (Augustus, MAKER, Funannotate, NCBI RefSeq) mixed in one analysis. **Mechanism:** Pipelines differ in handling of intron predictions, gene boundaries, and small-gene filtering. A gene predicted by Augustus but missed by MAKER appears as an apparent lineage-specific gene in the MAKER-annotated species, producing spurious "species-specific" orthogroups. **Symptom:** CAFE (gene family evolution) shows extreme expansions / contractions for species with different annotation pipelines; per-species "unique" gene count varies 5-10x between technically similar genomes. **Fix:** Re-annotate all genomes with one pipeline (e.g. BRAKER3 or Funannotate) before orthology. When re-annotation isn't possible, normalize via BUSCO/Compleasm completeness (filter OGs absent from species with < 95% BUSCO complete). Document annotation pipeline per species in methods. ### Ortholog conjecture violations **Trigger:** Transferring GO/KEGG functional annotation from 1:1 ortholog without testing functional divergence. **Mechanism:** Studies show orthologs are *weakly* more functionally similar than paralogs (Altenhoff 2012 PLoS Comp Biol 8:e1002514), though Stamboulian 2020 (Bioinformatics 36:i219) finds paralogs often comparable predictors; effect size is small and dependent on evolutionary distance. Subfunctionalization (Force 1999 Genetics 151:1531), neofunctionalization, and dosage subfunctionalization can rapidly differentiate 1:1 orthologs. **Symptom:** Transferred annotation contradicts species-specific experimental data; ortholog has fold-change different expression patterns; rapid evolution (dN/dS > 0.3) on one branch only. **Fix:** For high-confidence annotation transfer, require: (1) 1:1 orthology AND (2) low branch dN/dS (< 0.2 for both branches) AND (3) GO evidence code ISO with experimental support upstream. Tag annotation as "predicted from ortholog" rather than equating function. eggNOG ortholog-conjecture-aware mode provides confidence scoring. ### RBH symmetric-but-wrong errors **Trigger:** Reciprocal best hits (RBH) method on a gene that has been replaced by a paralog in one lineage. **Mechanism:** Lineage A retains the original ortholog; Lineage B lost the ortholog and replaced its function with a paralog (xenologous replacement). RBH between A and B identifies the paralog as the "ortholog" because it's the best hit. **Symptom:** dN/dS on RBH ortholog pair > 1 on the B branch (suggesting positive selection but actually paralog substitution); gene tree shows the B sequence is sister to other paralogs of the A gene, not to A. **Fix:** Use tree-based methods (OrthoFinder3, OMA HOGs) for organisms with extensive paralog history. RBH is appropriate for orthogroup pre-screening only. JustOrthologs is explicitly an RBH method, so it's affected by this; reserve for closely related species. ### Synteny ignored in WGD lineages **Trigger:** Plant, yeast, fish, or salmonid analysis where one or more ancient WGD events occurred. **Mechanism:** WGD doubles all genes; subsequent gene loss is biased toward certain functional categories (Dosage Balance Hypothesis; Birchler & Veitia 2007 Plant Cell 19:395). Sequence-only orthology cannot distinguish "ortholog" (single ancestral gene) from "homeolog" (paralog from WGD). In hexaploids (wheat, Brassica), this triples the confusion. **Symptom:** Multiple co-orthologs per species in orthogroups for WGD-affected lineages; rate variation across "co-orthologs" suggests one is the true ortholog and the others are recent duplicates. **Fix:** Use synteny-aware methods (GENESPACE -- Lovell 2022 eLife 11:e78526; ProteinOrtho `-synteny`; OrthoFinder3 + post-hoc synteny verification). Restrict 1:1 ortholog phylogenomics to genes in single-copy syntenic regions. For deep WGD ancestry (e.g. 2R vertebrate WGD), modern orthologs of post-2R paralogs (ohnologs) can no longer be unambiguously identified by sequence; use synteny + duplication-dating. ### MAFFT-only alignment in OrthoFinder MSA mode **Trigger:** `orthofinder -M msa -A mafft` default; very divergent orthogroups (deep taxonomy). **Mechanism:** MAFFT default settings choose accurate-but-fast iterative refinement. For sequences with >50% divergence, MAFFT's auto choice may downgrade to FFT-NS-2, producing alignments with > 30% poorly aligned columns that distort gene-tree inference and downstream HOG construction. **Symptom:** Per-OG MAFFT logs show "FFT-NS-2 selected"; gene trees have unstable rooting; bootstrap support < 60% for many internal branches. **Fix:** Force `mafft-linsi` or specify `-A mafft --thread -1 --localpair --maxiterate 1000`. For deeply divergent OGs, run PRANK or MUSCLE5 post-hoc on critical OGs and re-build gene trees with IQ-TREE. PREQUAL or HmmCleaner segment-filtering improves resulting trees (Di Franco 2019 BMC Evol Biol 19:21). ## Quantitative Thresholds | Quantity | Threshold | Source / Rationale | |----------|-----------|-------------------| | QfO benchmark adoption | 100% within-species ortholog pairs identified | Altenhoff 2024 NAR Genom Bioinform 6:lqae167; minimum competence threshold | | BUSCO/Compleasm completeness for inclusion | >= 90% complete (single + duplicated) before orthology run | Below this, expect inflated lineage-specific OGs | | OrthoFinder3 outgroup distance | >= one sister taxonomic level (sister phylum / class / order) | Emms & Kelly 2017 MBE 34:3267 STRIDE rooting | | Single-copy ortholog filter for phylogenomics | Present in >= 90% of species, exactly 1 copy each | Standard convention; below this, missing data biases tree inference | | MMseqs2 sensitivity for divergent homologs | `-s 7.5` for > 50% divergence | mmseqs2 documentation; default `-s 4.0` misses distant homologs | | ProteinOrtho `--conn` (connectivity) | >= 0.1 default; 0.2 stricter for less paralogy confusion | Lechner 2011 | | OMA HOG inclusion criterion | RBH + Smith-Waterman score and pairwise stability | OMA convention; sub-clade HOGs nested in supergroup HOGs | | Ortholog age threshold for functional transfer | divergence < 200 Myr OR dN/dS < 0.2 | Operational convention (ortholog conjecture weaker at higher divergence) | | Annotation pipeline normalization | All species annotated with same pipeline OR BUSCO-completeness within 5% | Avoid annotation-heterogeneity bias on CAFE | | eggNOG-mapper minimum score | bit score / e-value defaults; check `--seed_ortholog_score` | eggNOG-mapper docs | | TOGA intactness classes (loss_summ_data.tsv) | I (intact), PI (partial intact), UL (uncertain loss), L (lost), M (missing/assembly gap), PM (partial missing) | Kirilenko 2023 + TOGA repo | | TOGA orthology relationships (orthology_classification.tsv) | one2one, one2many, many2one, many2many, PG (paralogous projection) | Kirilenko 2023 | | SonicParanoid2 confidence | >= 0.9 high-confidence orthologs; 0.5-0.9 moderate | Cosentino 2024 docs | | Reasonable expected runtime (200 species, 16 cores) | OrthoFinder3 default ~8-24h; SonicParanoid2 ~2-6h; FastOMA ~3-8h | Hardware-dependent benchmarks | ## OrthoFinder3 Standard Workflow **Goal:** Produce HOG-based orthology with species tree and gene-duplication events for any clade. **Approach:** Prepare cleaned per-species proteomes (longest isoforms) -> include outgroup -> run with MSA + IQ-TREE option -> parse HOG output at appropriate taxonomic level. ```bash # Pre-clean proteomes: longest isoform per gene for f in raw_proteomes/*.faa; do python tools/primary_transcript.py $f > cleaned/$(basename $f) done # OrthoFinder v3 with MSA + IQ-TREE for tree-based HOG inference orthofinder \ -f cleaned/ \ -t 16 \ -a 4 \ -M msa \ -A mafft \ -T iqtree \ -S diamond_ultra_sens \ -y \ -o orthofinder_run # Output of interest (v3 layout): # orthofinder_run/Results_<date>/Phylogenetic_Hierarchical_Orthogroups/N0.tsv (root-level HOGs) # orthofinder_run/Results_<date>/Single_Copy_Orthologue_Sequences/ (single-copy MSA-ready) # orthofinder_run/Results_<date>/Species_Tree/SpeciesTree_rooted.txt # orthofinder_run/Results_<date>/Gene_Duplication_Events/ (per-branch dup counts) ``` ```python '''Parse OrthoFinder v3 HOG output for downstream analysis.''' import pandas as pd from pathlib import Path def load_hogs(results_dir, level='N0'): '''HOG file path differs from v2: Phylogenetic_Hierarchical_Orthogroups/{level}.tsv.''' p = Path(results_dir) / 'Phylogenetic_Hierarchical_Orthogroups' / f'{level}.tsv' df = pd.read_csv(p, sep='\t') return df.set_index('HOG') def _is_present(v): '''OrthoFinder v3 HOG cells are either NaN (older versions), '' (newer), or a comma-separated gene list.''' return not (pd.isna(v) or v == '' or str(v).strip() == '') def single_copy_hogs(hog_df, min_species_fraction=0.9): '''Return HOG IDs where every present species has exactly 1 gene and >= fraction of species present.''' sp_cols = [c for c in hog_df.columns if c not in ('OG', 'Gene Tree Parent Clade')] is_single = hog_df[sp_cols].apply( lambda r: all((not _is_present(v)) or ',' not in str(v) for v in r), axis=1 ) n_present = hog_df[sp_cols].apply(lambda r: sum(_is_present(v) for v in r), axis=1)