Multi-species orthology with OrthoFinder

Comparative genomics

scRNA-seq

HPC & computing

Find orthogroups and pairwise orthologs across several cnidarian proteomes in a single OrthoFinder run, ready for cross-species marker transfer and comparative genomics. One LSF job.

Published

May 25, 2026

Modified

October 16, 2026

Pipeline·Tools OrthoFinder 3 MAFFT FastTree DIAMOND

Overview

OrthoFinder groups the proteins of several species into orthogroups (gene families descended from one ancestral gene) and uses gene trees to call orthologs between every pair of species. In the single-cell work these tables are the bridge between species: published cell-type markers from one coral atlas are mapped onto genes in another species to score its clusters.

Run all species together, not pairwise

Coral genomes have many tandem duplications. Gene trees built from several species at once (plus an outgroup such as Nematostella vectensis) root correctly and separate orthologs from paralogs much better than separate two-species runs. All recent analyses use one combined run.

flowchart LR
  P[Proteome FASTA or URL<br/>per species] --> C[Clean + check<br/>login node]
  C --> O[OrthoFinder<br/>MSA gene trees<br/>one LSF job]
  O --> G[(Orthogroups.tsv<br/>gene counts<br/>single-copy OGs)]
  O --> R[(Pairwise ortholog<br/>tables)]
  R --> S[Pair summary:<br/>total and 1:1 pairs]

Inputs

Input	Notes
One protein FASTA per species	Local path or an `https://` URL (e.g. NCBI, reefgenomics, GitHub). Gzipped files are fine.
A short species code per species	Becomes the file name (`Pdam.fa`) and therefore the column name in `Orthogroups.tsv`. Letters, numbers and `_` only.
An outgroup	e.g. N. vectensis (`NV2g.20240221.protein.fa` from SimrBase)

Ideally use one protein per gene (longest isoform). OrthoFinder treats every isoform as a separate gene, which inflates paralogs.

Software

Tool	How it’s loaded
OrthoFinder (with MAFFT and FastTree for MSA gene trees)	Its own conda env, e.g. `conda create -n orthofinder_env -c bioconda orthofinder`

Keep orthofinder_env separate from other environments (see the functional annotation pipeline, which uses the same rule).

Run it

The script never needs editing; your species and paths go in a config file.

Make a config file:

bash orthofinder_multispecies.sh --init my_orthofinder.config

Fill it in (nano my_orthofinder.config). Replace every <fill in: ...> value:

RUN="<fill in: short run name, e.g. Pdam_Amil_Spis_Nvec>"
RUN_DIR="<fill in: /path/to/scratch/orthofinder/RUN_NAME>"
SPECIES=(   "<fill in: code 1, e.g. Pdam>"   "<fill in: code 2>"   "<fill in: code 3>" )
PROTEOMES=( "<fill in: /path/or/URL/to/proteins_1.fasta>" "<fill in: proteins_2>" "<fill in: proteins_3>" )
# ... plus LSF project, email; defaults for cores, memory and OrthoFinder options

Check the proteomes first (fast, on the login node):
```
PREP_ONLY=1 bash orthofinder_multispecies.sh my_orthofinder.config
```
This writes cleaned files to <RUN_DIR>/proteomes/<CODE>.fa and a proteome_summary.tsv with protein counts. It also stops if a download fails, a file has no sequences, or protein IDs are duplicated after shortening headers. Cleaning removes any text before the first >, removes blank lines, strips Windows line endings, and shortens each header to its first word.
Submit:
```
bash orthofinder_multispecies.sh my_orthofinder.config
bjobs <JOB ID>
```
Default resources: bigmem, 16 cores, rusage[mem=15000], 120 h walltime (the settings used for the four-species run below).

Full script

Download orthofinder_multispecies.sh

orthofinder_multispecies.sh

#!/bin/bash
# =============================================================================
# Multi-species orthology with OrthoFinder (LSF clusters, e.g. Pegasus)
# Pipeline page: https://kevinhwong1.github.io/KevinHWong_Notebook/pipelines/orthofinder-multispecies.html
#
# What it does:
#   1. Prepares one clean protein FASTA per species (download if a URL is given; drop any text before
#      the first header and blank lines; shorten headers to the first word; check for duplicate IDs).
#      Each file is named <SPECIES>.fa, so species codes become the column names in Orthogroups.tsv.
#   2. Submits ONE OrthoFinder job with all species together (MSA gene trees by default).
#   3. When it finishes, copies the key tables to <RUN_DIR>/summary/ and counts ortholog pairs.
#
# Usage:
#   bash orthofinder_multispecies.sh --init my_orthofinder.config   # 1. write a config template
#   nano my_orthofinder.config                                      # 2. replace every placeholder
#   bash orthofinder_multispecies.sh my_orthofinder.config          # 3. prepare FASTAs + submit
#   PREP_ONLY=1 bash orthofinder_multispecies.sh my_orthofinder.config  # only prepare/check FASTAs
#
# Nothing in this file needs editing; all run-specific settings live in the config file.
# =============================================================================
set -euo pipefail

usage() { sed -n '13,17p' "$0" | sed 's/^# \{0,1\}//'; }

write_template() {
cat > "$1" <<'TEMPLATE'
# Config for orthofinder_multispecies.sh
# Replace every "<fill in: ...>" value. The script refuses to run while any are left.
# This file is sourced by bash: no spaces around "=", and quote values.

# --- Run and cluster ----------------------------------------------------------
RUN="<fill in: short run name, e.g. Pdam_Amil_Spis_Nvec>"
RUN_DIR="<fill in: /path/to/scratch/orthofinder/RUN_NAME>"
LSF_PROJECT="<fill in: LSF project/allocation passed to bsub -P>"
EMAIL="<fill in: email for LSF job notifications>"
QUEUE="bigmem"
CORES=16
MEM_MB=15000
WALLTIME="120:00"

# --- Species: one entry per species, same order in both lists ----------------------
# SPECIES codes become the column names in Orthogroups.tsv (short, no spaces, e.g. Pdam).
# PROTEOMES can be a local path or an https:// URL (downloaded on the login node).
# Include an outgroup (e.g. Nematostella vectensis) to help OrthoFinder root gene trees.
SPECIES=(   "<fill in: code 1, e.g. Pdam>"   "<fill in: code 2>"   "<fill in: code 3>" )
PROTEOMES=( "<fill in: /path/or/URL/to/proteins_1.fasta>" "<fill in: proteins_2>" "<fill in: proteins_3>" )

# --- Software ---------------------------------------------------------------------
CONDA_SH="${HOME}/anaconda3/etc/profile.d/conda.sh"
ORTHOFINDER_ENV="orthofinder_env"     # OrthoFinder + MAFFT + FastTree; keep separate from other envs

# MSA gene trees (more accurate orthologs than the default DendroBLAST trees, but slower)
ORTHOFINDER_ARGS="-M msa -A mafft -T fasttree"
TEMPLATE
}

# ----------------------------------------------------------------------------- read config
if [[ "${1:-}" == "--init" ]]; then
  cfg="${2:-orthofinder.config}"
  if [[ -e "${cfg}" ]]; then echo "${cfg} already exists; not overwriting." >&2; exit 1; fi
  write_template "${cfg}"
  echo "Wrote ${cfg}. Fill in every <fill in: ...> value, then run:"
  echo "  bash $0 ${cfg}"
  exit 0
fi

CONFIG="${1:-}"
if [[ -z "${CONFIG}" || ! -f "${CONFIG}" ]]; then usage; exit 1; fi
if grep -nE '^[^#]*<fill in' "${CONFIG}"; then
  echo "" >&2; echo "Replace the placeholder values above in ${CONFIG} first." >&2; exit 1
fi
# shellcheck source=/dev/null
source "${CONFIG}"

# ----------------------------------------------------------------------------- checks
if [[ ${#SPECIES[@]} -ne ${#PROTEOMES[@]} ]]; then
  echo "SPECIES and PROTEOMES must have the same number of entries." >&2; exit 1
fi
if [[ ${#SPECIES[@]} -lt 2 ]]; then echo "Give at least 2 species." >&2; exit 1; fi
if [[ $(printf '%s\n' "${SPECIES[@]}" | sort | uniq -d | wc -l) -gt 0 ]]; then
  echo "SPECIES codes must be unique." >&2; exit 1
fi
for s in "${SPECIES[@]}"; do
  [[ "${s}" =~ ^[A-Za-z0-9_]+$ ]] || { echo "Species code '${s}': use only letters, numbers and _." >&2; exit 1; }
done

FASTA_DIR="${RUN_DIR}/proteomes"
RESULTS_PARENT="${RUN_DIR}/results"
RESULTS_DIR="${RESULTS_PARENT}/Results_${RUN}"
LOGDIR="${RUN_DIR}/logs"
SUMMARY_DIR="${RUN_DIR}/summary"
mkdir -p "${FASTA_DIR}" "${LOGDIR}"

if [[ -e "${RESULTS_DIR}" ]]; then
  echo "${RESULTS_DIR} already exists. Use a new RUN name (or move the old results)." >&2; exit 1
fi

# ----------------------------------------------------------------------------- 1. prepare FASTAs
echo "== Preparing proteomes in ${FASTA_DIR} =="
printf "species\tsequences\tsource\n" > "${RUN_DIR}/proteome_summary.tsv"
problems=0
for i in "${!SPECIES[@]}"; do
  sp="${SPECIES[$i]}"; src="${PROTEOMES[$i]}"; out="${FASTA_DIR}/${sp}.fa"
  tmp="${FASTA_DIR}/.${sp}.download"
  if [[ "${src}" =~ ^https?:// ]]; then
    echo "  ${sp}: downloading ${src}"
    curl -fsSL -o "${tmp}" "${src}" || { echo "  ${sp}: download failed" >&2; problems=$((problems+1)); continue; }
    in="${tmp}"
  else
    [[ -f "${src}" ]] || { echo "  ${sp}: not found: ${src}" >&2; problems=$((problems+1)); continue; }
    in="${src}"
  fi
  # unzip if needed; drop text before the first '>' and blank lines; keep the first word of each header
  if [[ "${in}" == *.gz ]] || gzip -t "${in}" 2>/dev/null; then cat_cmd="gzip -dc"; else cat_cmd="cat"; fi
  ${cat_cmd} "${in}" | tr -d '\r' | awk '
    /^>/ { seen=1; split(substr($0,2), a, /[ \t]/); print ">" a[1]; next }
    seen && NF { print }' > "${out}"
  rm -f "${tmp}"

  n=$(grep -c '^>' "${out}" || true)
  dups=$(grep '^>' "${out}" | sort | uniq -d | wc -l)
  if (( n == 0 )); then echo "  ${sp}: no sequences found!" >&2; problems=$((problems+1)); fi
  if (( dups > 0 )); then echo "  ${sp}: ${dups} duplicate protein IDs after shortening headers" >&2; problems=$((problems+1)); fi
  printf "%s\t%s\t%s\n" "${sp}" "${n}" "${src}" >> "${RUN_DIR}/proteome_summary.tsv"
  echo "  ${sp}: ${n} proteins"
done
if (( problems > 0 )); then echo "Fix the ${problems} problem(s) above before running OrthoFinder." >&2; exit 1; fi

echo "Check that each species has one protein per gene (longest isoform). If a proteome has many"
echo "isoforms per gene, filter it first; OrthoFinder treats every isoform as a separate gene."

if [[ "${PREP_ONLY:-0}" == 1 ]]; then echo "PREP_ONLY=1: stopping before submission."; exit 0; fi

# ----------------------------------------------------------------------------- 2. submit OrthoFinder
JOBFILE="${LOGDIR}/orthofinder_${RUN}.lsf"
cat > "${JOBFILE}" <<EOF
#!/bin/bash
#BSUB -J OF_${RUN}
#BSUB -q ${QUEUE}
#BSUB -P ${LSF_PROJECT}
#BSUB -n ${CORES}
#BSUB -W ${WALLTIME}
#BSUB -R "rusage[mem=${MEM_MB}]"
#BSUB -u ${EMAIL}
#BSUB -o ${LOGDIR}/orthofinder_${RUN}_%J.out
#BSUB -e ${LOGDIR}/orthofinder_${RUN}_%J.err
#BSUB -B
#BSUB -N
set -eo pipefail   # no -u: conda activate scripts reference unset variables

source ${CONDA_SH}
conda activate ${ORTHOFINDER_ENV}
orthofinder --help | head -2 || true

orthofinder -f ${FASTA_DIR} -o ${RESULTS_PARENT} -n ${RUN} -t ${CORES} -a ${CORES} ${ORTHOFINDER_ARGS}

# ---- summary tables
mkdir -p ${SUMMARY_DIR}
cp ${RESULTS_DIR}/Orthogroups/Orthogroups.tsv \\
   ${RESULTS_DIR}/Orthogroups/Orthogroups.GeneCount.tsv \\
   ${RESULTS_DIR}/Orthogroups/Orthogroups_SingleCopyOrthologues.txt \\
   ${RESULTS_DIR}/Comparative_Genomics_Statistics/Statistics_Overall.tsv \\
   ${SUMMARY_DIR}/

# ortholog pairs per species pair: total rows and strict 1:1 rows (one gene on each side)
printf "species_a\tspecies_b\tortholog_rows\tone_to_one\n" > ${SUMMARY_DIR}/ortholog_pair_summary.tsv
for f in ${RESULTS_DIR}/Orthologues/Orthologues_*/*__v__*.tsv; do
  pair=\$(basename "\$f" .tsv)
  awk -F'\t' -v a="\${pair%%__v__*}" -v b="\${pair##*__v__}" '
    NR>1 { rows++; if (\$2 !~ /,/ && \$3 !~ /,/) one++ }
    END  { printf "%s\t%s\t%d\t%d\n", a, b, rows, one }' "\$f" >> ${SUMMARY_DIR}/ortholog_pair_summary.tsv
done
echo "Done. Key tables in ${SUMMARY_DIR}/; full results in ${RESULTS_DIR}/"
EOF

JID=$(bsub < "${JOBFILE}" | grep -oP '(?<=Job <)\d+')
cat <<EOF

Submitted OrthoFinder job ${JID} for ${#SPECIES[@]} species (${SPECIES[*]}).
Monitor with:     bjobs ${JID}      (or bpeek ${JID})
Job script:       ${JOBFILE}
Results will be:  ${RESULTS_DIR}/
Key tables:       ${SUMMARY_DIR}/  (Orthogroups.tsv, ortholog_pair_summary.tsv, ...)
EOF

Outputs

In <RUN_DIR>/summary/ (copied there when the job finishes):

File	What it is
`Orthogroups.tsv`	Every orthogroup, with that species’ genes in each column
`Orthogroups.GeneCount.tsv`	Genes per species per orthogroup (gene-family expansions)
`Orthogroups_SingleCopyOrthologues.txt`	Orthogroups with exactly one gene in every species
`Statistics_Overall.tsv`	% genes assigned, number of orthogroups, G50/O50
`ortholog_pair_summary.tsv`	For each species pair: ortholog rows and strict 1:1 pairs

Pairwise ortholog tables for each species are in <RUN_DIR>/results/Results_<RUN>/Orthologues/Orthologues_<CODE>/<CODE>__v__<OTHER>.tsv.

Example results

Four species for P. damicornis cell-type annotation (P. damicornis, A. millepora, S. pistillata, N. vectensis; May 2026): 103,775 genes (89.3%) assigned to 17,897 orthogroups; 9,303 orthogroups contain all four species, and 4,546 of those are single-copy in every species.

Pair	Ortholog rows	1:1 pairs
P. damicornis – A. millepora	11,797	8,657
P. damicornis – S. pistillata	16,788	12,937
P. damicornis – N. vectensis	11,410	6,683

Ten species for the cnidarian single-cell atlas (A. cervicornis, A. millepora, A. muricata, G. fascicularis, N. vectensis, O. patagonica, O. faveolata, P. damicornis, S. pistillata, Xenia sp.): 19,320 orthogroups, 1,146 single-copy across all ten.

Using the tables for marker transfer

Gene IDs rarely match exactly between a published atlas, its proteome and your own Cell Ranger features. Before scoring clusters with transferred markers:

Normalize IDs on each side (e.g. Pdam: CellBender pdam_00025168 → pdam-00025168; proteome/OrthoFinder pdam_00000427-RA; Gfas: strip the trailing .m1).
Check the match rate empirically: count how many marker-table IDs are actually found in Orthogroups.tsv after normalization. In the atlas this was 97.8% (A. millepora), 98.7% (S. pistillata) and 96.6% (O. patagonica). For A. muricata the marker IDs (EVM models) and the proteome IDs (RefSeq) had no crosswalk, so that reference was left out of automated scoring.
Prefer 1:1 orthologs for module scores; one-to-many groups are better for manual follow-up.

Known issues

Messy FASTA files from GitHub and other sources

The A. millepora and S. pistillata proteomes downloaded from GitHub had a trailing blank line, and one had text before the first header. Long headers with spaces and pipes also cause trouble downstream. The prep step handles all of these and stops on duplicate IDs.

Results folder already exists

OrthoFinder won’t write into an existing results folder. The script names results Results_<RUN> and stops before submitting if that folder exists. Use a new RUN name.

Default OrthoFinder vs MSA gene trees

The default (DendroBLAST) is much faster, and is fine for a first look. MSA gene trees (-M msa -A mafft -T fasttree) are slower but give better ortholog calls in gene families with many duplications, so they’re the default here. To use the fast mode, set ORTHOFINDER_ARGS="" in the config.

How it connects

Coral proteome functional annotation: also runs OrthoFinder per species (vs human, mouse and N. vectensis) to get best human orthologs.
scRNA-seq QC and clustering: clusters from here are annotated using markers transferred through these orthologs.
Earlier notebook entries: Troubleshooting OrthoFinder, OrthoFinder between coral species for cell type annotation.

References

Emms DM & Kelly S 2019. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biology 20:238.
Levy S et al. 2025. Nature 648:368–376 (source of the A. millepora, S. pistillata and O. patagonica proteomes and atlases).

Changelog

Date	Change
2026-10-05	Config-driven script: `--init` template, FASTA cleaning and checks, `Results_<RUN>` naming, pair summary
2026-08	Ten-species combined run replaces separate pairwise tables for the atlas
2026-05	Four-species run with MSA gene trees for P. damicornis annotation
2025-11	Pairwise runs (Acer–Amil, Pdam–Spis, Acer–Pdam–Nvec)
2023-12	First OrthoFinder tests