Bio Comparative Genomics Gene Family Evolution

Version Compatibility

Reference examples tested with: CAFE5 5.1.0+ (Mendes et al 2020 Bioinformatics 36(22-23):5516-5518), Count 11.0319+ (Csurös 2010 Bioinformatics 26:1910), BadiRate 1.35+ (Librado 2012 Bioinformatics 28:279), DupliPHY-Family (Ames et al 2012), CAFExp (legacy CAFE 4.2 -- DEPRECATED; use CAFE5), OrthoFinder 3.0+ for HOG input, R 4.4+, mclust 6.1+, phytools 2.3+, ETE4 4.1.0+ for tree manipulation. ALE/GeneRax/AleRax in companion skill [[gene-tree-species-tree-reconciliation]].

Before using code patterns, verify installed versions match. If versions differ:

• CLI: cafe5 --help; Count.exe (Java); badirate --help
• R: packageVersion('phytools')
• Python: pip show ete4

If code throws CAFE5: lambda did not converge, Count negative branch length, BadiRate gamma not initialized, the most common causes are: (1) annotation heterogeneity inflating family sizes, (2) saturated families (CAFE5 needs reasonable rate variation), (3) negative branch lengths in input tree (Count requires ultrametric). Pre-process: filter OG matrix to families present in >= 50% of species; resolve polytomies; ultrametricize tree.

Gene Family Evolution

"Which gene families expanded or contracted in which lineages?" -> Birth-death models on phylogeny (Hahn 2005; Csurös 2010) treat each orthogroup's per-species count as evolving under a stochastic birth-death process; lineage-specific rate shifts are detected as departures from a global rate. Annotation heterogeneity is the single largest confounder: different annotation pipelines predict different numbers of genes per family, producing apparent lineage-specific expansions that are artifacts of annotation choice (Tonkin-Hill 2020 demonstrated this for bacterial pangenomes). Consistent annotation + BUSCO/Compleasm completeness filtering are mandatory before any birth-death model interpretation. CAFE5 (Mendes et al 2020 Bioinformatics 36:5516) replaces older CAFE versions with gamma-distributed rate categories for more biologically realistic modeling.

• CLI: cafe5 -i orthogroup_counts.tsv -t species_tree.nwk -p -- main CAFE5 workflow
• CLI: cafe5 -e -- error-aware mode for annotation uncertainty
• CLI: Count -- Java GUI / CLI for parsimony + likelihood ASR
• CLI: badirate -t tree.nwk -d counts.tsv -- likelihood birth-death + branch parsimony

Algorithmic Taxonomy

Methodology evolves; CAFE5 is the modern standard; verify the current CAFE5 manual (hahnlab/CAFE5) and Hahn lab papers before locking on a single approach.

Decision Tree by Experimental Scenario

Per-Tool Failure Modes

Annotation heterogeneity inflating expansions / contractions

Trigger: Running CAFE5 on counts from genomes annotated by different pipelines (Augustus, MAKER, BRAKER, NCBI).

Mechanism: Different annotation tools predict different numbers of genes per family; the same biological gene family may be annotated with 5 genes in BRAKER and 8 in MAKER. CAFE5 sees the difference as an "expansion" in the MAKER-annotated species (Tonkin-Hill 2020 documented this for bacterial pangenomes; same principle in eukaryotes).

Symptom: "Most expanded families" cluster in species annotated by a single pipeline (often more permissive tool); per-species "expansion rate" correlates with annotation pipeline rather than biology.

Fix: Re-annotate all genomes with a single pipeline (currently BRAKER3 or Funannotate for eukaryotes; Bakta for bacteria) before CAFE5. Alternatively use CAFE5-error mode with explicit error rates per species. Document annotation pipeline + version per species.

Assembly fragmentation creating false contractions

Trigger: Including draft assemblies with low N50 in CAFE5 analysis.

Mechanism: Fragmented assemblies miss genes; the same family appears with fewer genes than expected. CAFE5 reports this as a contraction in the affected species.

Symptom: "Contracted" families in fragmented assemblies; correlation between BUSCO completeness and CAFE5 "contraction" rate; per-species missing-gene count varies 5-10x.

Fix: Require >= 90% BUSCO/Compleasm completeness for inclusion. Exclude species with > 5% lower BUSCO than median. Document N50 + BUSCO per assembly. For unavoidably fragmented assemblies, exclude from CAFE5 or use CAFE5-error with empirical error estimates.

CAFE5 lambda non-convergence

Trigger: CAFE5 reports "lambda did not converge"; lambda jumps between values across runs.

Mechanism: Insufficient data (< 100 families); strong rate heterogeneity not captured; or input tree non-ultrametric.

Symptom: lambda estimate unstable; AIC of model selection variable.

Fix: Require >= 100 orthogroups (preferably > 1000) in input. Ensure tree is ultrametric (ape::chronos() or treePL); CAFE5 expects time-scaled tree. Use gamma rate categories (CAFE5 default -k 4) for rate heterogeneity. If still non-convergent, restrict to single-copy or small families; check tree branch lengths for negative values.

Gamma rate-category misinterpretation

Trigger: Reporting "gamma rate categories" as biological gene-family clusters.

Mechanism: CAFE5 gamma categories are statistical buckets representing rate heterogeneity across families; they're not "fast-evolving family clusters" with biological meaning per se.

Symptom: Confusion in interpretation; "category-1 families are special."

Fix: Treat gamma categories as a statistical device. Report per-family lambda (estimated under per-family rate model) or per-clade lambda. Functional interpretation comes from family-level statistical tests, not category assignment.

Multiple-testing across many families

Trigger: Reporting "significant" expansions / contractions without correction.

Mechanism: With ~10,000 families tested, ~500 will be significant at p=0.05 under H0. Without correction, false-discovery rate is high.

Symptom: Long lists of "expanded" families; functional enrichment dominated by chance hits.

Fix: Apply FDR (Benjamini-Hochberg) across families; or restrict to a priori hypothesized families. CAFE5 reports per-family p-values; FDR-correct downstream.

Tree non-ultrametric / negative branches

Trigger: Using ML tree directly (substitutions per site) as input to CAFE5.

Mechanism: CAFE5 expects ultrametric (time-calibrated) tree; substitution-based branch lengths are not time-calibrated and may even have negative branches after rate variation.

Symptom: CAFE5 errors with "negative branch length"; or produces unreliable lambda estimates.

Fix: Time-calibrate the tree: use ape::chronos(), treePL (Smith & O'Meara 2012 Bioinformatics 28:2689), or LSD2 (To et al 2016 Syst Biol 65:82; gascuel-lab/LSD2) for fast NP-like calibration. Or use existing time-calibrated tree from TimeTree database.

Outlier-family-driven lambda estimate

Trigger: Including extreme-size families (e.g., NLR resistance gene clusters in plants with 500+ members) in CAFE5.

Mechanism: Birth-death model's likelihood is dominated by large families; one or two outlier families can overwhelm the global lambda estimate.

Symptom: Excluding the top-5 largest families changes lambda by > 30%; lambda confidence intervals huge.

Fix: Robust analysis: report lambda with and without outliers; consider per-family lambda for largest families. Functional annotation of outliers reveals if they're biologically expected expansions (rapidly evolving gene families).

Convergent gene-family-rate shifts not captured

Trigger: Phylogenomic question is whether multiple independent lineages show similar gene-family-rate shifts.

Mechanism: CAFE5 doesn't natively test for convergent rate shifts; per-clade lambda is the closest analog.

Symptom: Manual inspection of expansions across independent lineages shows pattern; CAFE5 doesn't formalize it.

Fix: Combine CAFE5 per-family lambdas with RERconverge (Redlich et al 2024 MBE 41:msae210) for trait-correlated rate shifts; use CSUBST (Fukushima 2023 Nat Eco Evo 7:155) for convergent substitution patterns.

CAFE5 vs ALE for HGT-affected family

Trigger: Applying CAFE5 to bacterial families where transfer is common.

Mechanism: CAFE5 models birth-death (D, L) only; ignores transfer (T). For HGT-affected families, the dynamics include T events that CAFE5 cannot capture.

Symptom: CAFE5 lambda for bacterial families is implausibly high; family counts don't fit birth-death model.

Fix: Use ALE / GeneRax / AleRax for HGT-affected families ([[gene-tree-species-tree-reconciliation]]); CAFE5 is appropriate for vertical-inheritance-dominated families. Combine: CAFE5 for the bulk; ALE for HGT-confirmed families.