Version Compatibility

Reference examples tested with: BioPython 1.83+, BUSCO 5.5+, NCBI BLAST+ 2.15+, OrthoFinder 2.5+, pandas 2.2+

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

• Python:

  pip show <package>

  then

  help(module.function)

  to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Ortholog Inference

"Find orthologs across my species" → Identify orthologous gene groups, paralogs, and co-orthologs across multiple species using sequence similarity clustering and gene tree reconciliation.

• CLI:

  orthofinder -f proteomes/

  for all-vs-all orthogroup inference

Method Selection

Method	Approach	Best For	Tradeoff
OrthoFinder	Tree-based (gene tree reconciliation with species tree)	Accuracy, evolutionary analysis, gene duplication events	Slower, needs sufficient species
ProteinOrtho	Graph-based (reciprocal best hits + connectivity)	Speed, many genomes, quick surveys	Less accurate for complex gene families
OMA/FastOMA	Graph-based (strict pairwise, hierarchical groups)	Precision-critical applications, large-scale (1000+ genomes)	Lowest recall (misses distant orthologs)
SonicParanoid2	Graph-based (ML predictor + protein language model)	Fast + accurate graph-based	Newer, less community testing

Tree-based methods (OrthoFinder) build gene trees and reconcile with the species tree to distinguish speciation (orthology) from duplication (paralogy). More accurate but computationally expensive.

Graph-based methods (ProteinOrtho, OMA, SonicParanoid) use sequence similarity with clustering. Faster but can confuse paralogs with orthologs when evolutionary rates vary.

Default recommendation: OrthoFinder for most analyses. ProteinOrtho for quick surveys or 50+ genomes. OMA/FastOMA when precision is paramount.

Input Quality

Annotation quality directly affects orthology inference. Heterogeneous annotations across species spuriously inflate lineage-specific gene counts, creating false gene family expansions/contractions in downstream CAFE analysis.

• Use consistent annotation pipelines across species when possible
• Verify proteome completeness with BUSCO/Compleasm before running orthology
• Remove isoforms (keep longest per gene) to avoid inflating copy numbers
• Incomplete gene models produce truncated proteins that split true orthogroups

Orthology Subtypes

• One-to-one orthologs: single gene in each species, ideal for phylogenomics
• One-to-many / many-to-many: lineage-specific duplications after speciation
• In-paralogs: paralogs from duplication AFTER the speciation event of reference
• Out-paralogs: paralogs from duplication BEFORE the speciation event
• Co-orthologs: in-paralogous genes collectively orthologous to a gene in the outgroup

OrthoFinder Workflow

Goal: Infer orthologous gene groups across multiple species from their proteomes.

Approach: Run OrthoFinder on a directory of per-species FASTA files to perform all-vs-all DIAMOND search, gene/species tree inference, and ortholog/paralog classification, then parse the resulting orthogroups and classify by copy number pattern.

'''Ortholog inference with OrthoFinder'''

import subprocess
import pandas as pd
import os


def run_orthofinder(proteome_dir, output_dir=None, threads=4):
    '''Run OrthoFinder on directory of proteomes

    Input: Directory with one FASTA file per species
    File naming: Species name derived from filename

    OrthoFinder pipeline:
    1. All-vs-all DIAMOND/BLAST search
    2. Gene tree inference per orthogroup
    3. Species tree inference (STAG/STRIDE)
    4. Gene tree rooting and reconciliation
    5. Ortholog/paralog classification via DLC model

    Key options:
    -M msa: Use MSA-based gene trees (more accurate, slower; recommended for <20 species)
    -M dendroblast: Distance-based trees (default, faster; sufficient for >20 species)
    -S diamond: Fast search (default)
    -S blast: More sensitive (use for divergent species or small proteomes)
    '''
    cmd = f'orthofinder -f {proteome_dir} -t {threads}'

    if output_dir:
        cmd += f' -o {output_dir}'

    # Add -M msa for MSA-based gene trees (more accurate for evolutionary analysis)
    result = subprocess.run(cmd, shell=True, capture_output=True, text=True)

    # Output location
    if output_dir:
        results_dir = output_dir
    else:
        # OrthoFinder creates Results_MonDD in proteome_dir
        results_dir = None
        for d in os.listdir(proteome_dir):
            if d.startswith('OrthoFinder/Results_'):
                results_dir = os.path.join(proteome_dir, d)
                break

    return results_dir


def parse_orthogroups(orthogroups_file):
    '''Parse OrthoFinder Orthogroups.tsv

    Columns: Orthogroup, Species1, Species2, ...
    Values: Gene IDs (comma-separated if multiple)

    Orthogroup types:
    - Single-copy: One gene per species (ideal for phylogenomics)
    - Multi-copy: Duplications in some lineages
    - Species-specific: Genes unique to one species
    '''
    df = pd.read_csv(orthogroups_file, sep='\t')
    df = df.set_index('Orthogroup')

    orthogroups = {}
    for og_id, row in df.iterrows():
        genes = {}
        for species in df.columns:
            cell = row[species]
            if pd.notna(cell) and cell:
                genes[species] = cell.split(', ')
            else:
                genes[species] = []
        orthogroups[og_id] = genes

    return orthogroups


def classify_orthogroups(orthogroups, species_list):
    '''Classify orthogroups by copy number pattern

    Categories:
    - single_copy: Exactly one gene per species (best for phylogenomics)
    - universal: Present in all species (possibly multicopy)
    - partial: Missing from some species
    - species_specific: Only in one species
    '''
    classification = {
        'single_copy': [],
        'universal': [],
        'partial': [],
        'species_specific': []
    }

    for og_id, genes in orthogroups.items():
        present_in = [sp for sp in species_list if genes.get(sp)]
        copy_counts = [len(genes.get(sp, [])) for sp in species_list]

        if len(present_in) == 1:
            classification['species_specific'].append(og_id)
        elif len(present_in) == len(species_list):
            if all(c == 1 for c in copy_counts):
                classification['single_copy'].append(og_id)
            else:
                classification['universal'].append(og_id)
        else:
            classification['partial'].append(og_id)

    return classification


def get_single_copy_orthologs(orthogroups_file):
    '''Extract single-copy orthologs for phylogenomics

    Single-copy orthologs are ideal because:
    - Clear 1:1 relationships
    - No paralogy complications
    - Suitable for concatenated alignments
    '''
    df = pd.read_csv(orthogroups_file, sep='\t')
    df = df.set_index('Orthogroup')

    single_copy = []
    for og_id, row in df.iterrows():
        is_single = True
        for species in df.columns:
            cell = row[species]
            if pd.isna(cell) or cell == '':
                is_single = False
                break
            if ',' in str(cell):
                is_single = False
                break
        if is_single:
            single_copy.append(og_id)

    return df.loc[single_copy]

Gene Trees and Reconciliation

def parse_gene_trees(gene_trees_dir):
    '''Load gene trees from OrthoFinder

    Gene trees show evolutionary history within orthogroups
    Duplication/loss events inferred by species tree reconciliation
    '''
    from Bio import Phylo
    import glob

    trees = {}
    for tree_file in glob.glob(f'{gene_trees_dir}/*.txt'):
        og_id = os.path.basename(tree_file).replace('_tree.txt', '')
        trees[og_id] = Phylo.read(tree_file, 'newick')

    return trees


def identify_paralogs(orthogroup, species):
    '''Identify in-paralogs within an orthogroup

    In-paralogs: Duplications AFTER speciation (within one lineage)
    Out-paralogs: Duplications BEFORE speciation (separate orthogroups)
    Multiple genes from same species in an orthogroup = in-paralogs

    Distinguishing in- vs out-paralogs requires the species tree context
    and depends on which speciation event is being considered.
    OrthoFinder resolves this via gene tree reconciliation.
    '''
    genes = orthogroup.get(species, [])
    if len(genes) > 1:
        return {
            'species': species,
            'paralogs': genes,
            'count': len(genes)
        }
    return None


def find_co_orthologs(orthogroups, gene_id, species):
    '''Find co-orthologs of a gene

    Co-orthologs: Multiple genes in one species that are
    all orthologous to a single gene in another species

    Result of gene duplication after speciation
    '''
    for og_id, genes in orthogroups.items():
        if gene_id in genes.get(species, []):
            co_orthologs = {}
            for sp, sp_genes in genes.items():
                if sp != species and sp_genes:
                    co_orthologs[sp] = sp_genes
            return {'orthogroup': og_id, 'co_orthologs': co_orthologs}

    return None

ProteinOrtho Alternative

Goal: Detect orthologs using ProteinOrtho as a faster alternative for many-genome comparisons.

Approach: Run ProteinOrtho with DIAMOND backend on multiple proteome FASTA files and parse the output table for orthologous groups with connectivity scores.

def run_proteinortho(proteome_files, output_prefix, threads=4):
    '''Run ProteinOrtho for ortholog detection

    Faster than OrthoFinder for many genomes
    Uses synteny information if available

    -p=blastp+: Use DIAMOND (faster)
    -conn: Connectivity threshold (default 0.1)
    '''
    files_str = ' '.join(proteome_files)
    cmd = f'proteinortho -cpus={threads} -project={output_prefix} {files_str}'

    subprocess.run(cmd, shell=True)

    return f'{output_prefix}.proteinortho.tsv'


def parse_proteinortho(ortho_file):
    '''Parse ProteinOrtho output

    Columns: # Species, Genes, Alg.-Conn., Species1, Species2, ...
    '''
    df = pd.read_csv(ortho_file, sep='\t')

    orthogroups = {}
    for i, row in df.iterrows():
        og_id = f'OG{i:06d}'
        n_species = row['# Species']
        conn = row['Alg.-Conn.']

        genes = {}
        for col in df.columns[3:]:
            val = row[col]
            if pd.notna(val) and val != '*':
                genes[col] = val.split(',')
            else:
                genes[col] = []

        orthogroups[og_id] = {
            'genes': genes,
            'n_species': n_species,
            'connectivity': conn
        }

    return orthogroups

Functional Annotation Transfer

def transfer_annotation(query_gene, orthologs, annotation_db):
    '''Transfer functional annotation via orthology

    Confidence hierarchy:
    - One-to-one orthologs: Highest confidence; direct functional equivalence
    - Co-orthologs: Transfer to all, but note potential sub/neofunctionalization
    - In-paralogs (recent duplicates): Transfer with caution; function may have diverged
    - Distant orthologs (dS > 2): Lowest confidence; verify with domain conservation

    GO evidence codes:
    - ISO: Inferred from Sequence Orthology (recommended for 1:1 orthologs)
    - IBA: Inferred from Biological Aspect of Ancestor (phylogenetic propagation)
    - IEA: Inferred from Electronic Annotation (automated, lower confidence)

    Synteny context (see synteny-analysis) increases transfer confidence
    for genes in conserved genomic neighborhoods.
    '''
    annotations = []

    for species, genes in orthologs.items():
        for gene in genes:
            if gene in annotation_db:
                ann = annotation_db[gene]
                annotations.append({
                    'source_gene': gene,
                    'source_species': species,
                    'annotation': ann,
                    'evidence': 'ISO'  # Sequence orthology
                })

    return annotations

Completeness Assessment

Before orthology analysis, verify proteome completeness with BUSCO or Compleasm:

# BUSCO: standard benchmark against OrthoDB single-copy orthologs
busco -i proteome.fasta -m proteins -l <lineage> -o busco_out

# Compleasm: 14x faster alternative using miniprot
compleasm run -a genome.fasta -l <lineage> -o compleasm_out

BUSCO categories: Complete (single-copy + duplicated), Fragmented, Missing. Expect >90% complete for well-assembled genomes. High duplication rates may indicate assembly collapse or recent WGD. Choose the most specific available lineage for the clade being compared.

Related Skills

• comparative-genomics/synteny-analysis - Synteny-based ortholog verification and context
• comparative-genomics/positive-selection - Selection analysis on ortholog alignments
• phylogenetics/modern-tree-inference - Build species trees from single-copy orthologs
• alignment/pairwise-alignment - Align orthogroup sequences
• genome-annotation/annotation-transfer - Transfer annotations via orthology