1. 为什么我建议你用 OrthoFinder CAFE5 这条流程做基因家族收缩扩张分析最常被问到的就是“到底用哪套软件”。市面上的方案确实不少早年间有人用 OrthoMCL 做直系同源聚类再用 CAFEComputational Analysis of gene Family Evolution算家族变化后来有人改用 ProteinOrtho也有一部分人直接用 TreeFam 或者 Ensembl Compara 的现成结果。但如果你要我在 2025 年给一个“从零开始、结果可靠、可复现性高”的组合我会毫不犹豫推荐OrthoFinder CAFE5。OrthoFinder 的核心优势在于它不依赖预置的物种树而是通过基因树的综合比对来推断直系同源组orthogroup。相比 OrthoMCL 那种“先做 all-vs-all BLAST再用 MCL 聚类”的老路子OrthoFinder 对旁系同源paralog的处理更细致尤其是在近期发生了全基因组复制WGD的物种里效果会明显好很多。它输出的Orthogroups.tsv直接就是“基因家族 × 物种”的计数矩阵而 CAFE5 需要的输入恰好就是这个东西两个软件衔接得特别顺。CAFE5 的前身是 CAFE 3.x 和 CAFE 4.x核心算法延续了 Hahn 实验室的随机出生-死亡模型stochastic birth-death model通过最大似然估计来推断每个基因家族在系统发育树内部节点上的祖先大小再计算每个分支上的扩张expansion和收缩contraction显著性。CAFE5 相比旧版最大的改进是能更好地处理大规模数据、支持多线程并行而且对输入矩阵的格式容错度更高不再像 CAFE 4 那样动不动就抛出难懂的报错。这套流程也特别适合“非模式物种”的研究者。我见过很多人拿到的基因组注释质量一般基因结构预测出来的蛋白序列不完整但 OrthoFinder 对序列质量的容忍度比想象中高。只要你在过滤和格式清洗阶段稍微上点心后续分析依然能跑出稳定结果。简单说这套流程适合以下人群刚拿到一个或多个新基因组想快速了解基因家族在各物种间的扩张收缩情况做比较基因组学需要用基因家族变化支撑“某物种适应了某种环境”这类结论想复现文献里的分析结果但被各种旧版软件依赖问题折磨到崩溃的人。这篇文章就是冲着一个目标去的让你跟着操作从原始蛋白序列文件一路跑到 CAFE5 出结果中间不迷路。2. 核心思路拆解OrthoFinder 和 CAFE5 是如何接力的2.1 OrthoFinder 到底帮你做了什么很多教程会把 OrthoFinder 形容成“BLAST 聚类工具”这个说法不太准确。OrthoFinder 的工作流程大致是对输入的所有蛋白序列做 all-vs-all 比对默认用 Diamond速度比 BLASTP 快几个数量级基于比对结果通过 MCL 聚类形成初始的 orthogroup对每个 orthogroup 构建基因树利用基因树和物种树的关系识别并校正旁系同源区分直系同源和旁系同源输出最终的直系同源组矩阵。这个设计思路非常聪明它没有把“所有序列直接聚类到一个组”当成终点而是进一步用系统发育信息来校验每个组内部的关系。举个例子一个基因家族在物种 A 里有两个拷贝由串联复制产生在物种 B 里只有一个拷贝。MCL 可能会把 A 的两个拷贝和 B 的一个拷贝聚到同一个 orthogroup但 OrthoFinder 会在基因树里识别出 A 的两个拷贝是旁系同源然后把它们正确标注出来同时不影响整个 orthogroup 的计数。这对后续 CAFE5 的意义很大因为 CAFE5 需要的是“每个物种在该家族里的基因总数”如果你把旁系同源和直系同源混在一起数家族大小的波动会被严重扭曲。2.2 CAFE5 的数学原理用大白话讲清楚CAFE5 基于的是出生-死亡模型。你可以把每个基因家族想象成一个种群基因复制就是“出生”基因丢失就是“死亡”。这个模型假设在整个进化过程中每个谱系上基因家族的扩张和收缩速率是恒定的或者按你设定的多速率模型变化然后通过系统发育树的拓扑结构、分支长度和每个物种当前的基因家族大小反向推断出祖先节点上的家族大小。这里有一个非常核心的指标family-wide P-value。CAFE5 会针对每个基因家族在给定模型下计算“观察到的家族大小分布”出现的概率。如果概率很小说明这个家族的演化速率显著偏离了背景速率即这个基因家族在某个分支上发生了显著的扩张或收缩。要注意CAFE5 的 P-value 是针对整个家族的不是针对某一条分支的。如果你想知道“哪一条分支发生了显著变化”需要再加上-t参数进行分支特异性检验或者事后处理。很多初学者第一次跑完看到significant列里全是N就以为结果没意义其实只是因为你没有用对参数。2.3 为什么这两者的组合如此顺滑OrthoFinder 生成的Orthogroups.tsv是一个制表符分隔的矩阵行是 orthogroup ID列是物种单元格里是该物种在这个家族里的基因数量比如Orthogroup SpeciesA SpeciesB SpeciesC OG0000000 3 1 0 OG0000001 0 5 2CAFE5 需要的输入格式本质上也是“家族 × 物种”的计数矩阵唯一的区别是要把第一列换成家族 ID并额外加一行物种 ID 和一行树文件对应的物种名。所以从 OrthoFinder 到 CAFE5中间只差一个awk脚本转换的工作。这种“无缝衔接”是我强烈推荐这条组合的核心原因——不需要自己写一堆脚本来拼数据出错概率大大降低。3. 环境准备与软件安装含依赖坑点说明3.1 硬件与操作系统建议先说结论一套 8 核 16 线程、32 GB 内存的 Linux 服务器是起步配置。如果你只是拿 5 个物种、每个物种 3 万条蛋白来做测试16 GB 内存也能跑但 OrthoFinder 在构建基因树时很容易把内存占满。CAFE5 虽然轻量但如果你同时分析上万个基因家族内存占用也会线性增长。操作系统方面Ubuntu 20.04 / 22.04 LTS 或者 CentOS 7/8 都行。OrthoFinder 和 CAFE5 都提供了预编译的二进制包开箱即用。3.2 安装 OrthoFinderOrthoFinder 的安装非常无脑直接去 GitHub Releases 下载 Linux 版解压即可wget https://github.com/davidemms/OrthoFinder/releases/download/2.5.5/OrthoFinder_2.5.5.tar.gz tar -zxvf OrthoFinder_2.5.5.tar.gz cd OrthoFinder_2.5.5 ./orthofinder --help注意OrthoFinder 依赖 Diamond 和 MCL。新版的预编译包已经内置了 Diamond但 MCL 需要自己装sudo apt-get install mcl # 或者 yum install mcl如果你用的是新版 OrthoFinder它还需要fastme来做基因树。这个工具一般也自动处理了但建议装一下binutils以确保fastme能正常运行。安装完用orthofinder -h测一下能正常弹出帮助信息就没问题。3.3 安装 CAFE5CAFE5 的发布页在 GitHub 的hahnlab/CAFE5仓库下载 release 包后解压编译wget https://github.com/hahnlab/CAFE5/archive/refs/tags/v5.1.0.tar.gz tar -zxvf v5.1.0.tar.gz cd CAFE5-5.1.0 make如果编译报错多半是缺少依赖库。CAFE5 依赖libsqlite3和zlibsudo apt-get install libsqlite3-dev zlib1g-dev编译完成后cafe5的可执行文件会生成在bin/目录下。为了后续方便建议把它加到环境变量echo export PATH/path/to/CAFE5-5.1.0/bin:$PATH ~/.bashrc source ~/.bashrc3.4 安装 MAFFT 和 FastTree强烈建议虽然 OrthoFinder 运行时不是必须的外部工具但如果你要进行下游分析比如构建基因树做热图、或者用单个基因家族序列跑树MAFFT 和 FastTree 能大幅提升效率。这两者安装也很简单sudo apt-get install mafft fasttree实测下来MAFFT 的--auto模式在多序列比对方面完爆 ClustalW我后来所有涉及比对的地方都统一用 MAFFT不再回头看旧工具。4. 输入数据准备格式、命名与过滤技巧4.1 蛋白序列文件的格式要求OrthoFinder 要求输入的每个物种的蛋白序列保存在一个独立的 FASTA 文件中文件扩展名是.fa、.fasta或.faa。把这些文件统一放进一个文件夹里比如proteomes/。有两点需要特别留意序列头FASTA header不能有空格。OrthoFinder 默认以第一个空格作为 ID 分隔符如果序列头写成gene1 description xxx那它只会保留gene1作为 ID描述部分会被忽略。这个其实问题不大但如果你后续要追踪具体基因最好保持 ID 干净。不要在序列 ID 里用冒号和竖线。Ensembl 的蛋白序列经常是ENSP00000384948.1 pep:known这种格式里面那个冒号在后续处理中会引发莫名其妙的问题。建议直接用sed把序列头清洗成只有基因 IDfor f in proteomes/*.faa; do sed -i s/ .*// $f done4.2 序列过滤低质量序列会拖垮整体结果不是所有预测出来的蛋白序列都适合直接输入。我在实操中一定会做两步过滤去掉长度过短的序列。比如长度小于 30 个氨基酸的“蛋白质”大概率是假基因或基因组组装垃圾保留它们只会增加 BLAST 噪声。用seqkit一次性搞定seqkit seq -m 30 -g proteomes/speciesA.faa -o proteomes/speciesA.clean.faa用seqkit stat看每个文件里序列条数筛选掉极端异常的样本。比如某个物种注释出了 8 万个基因而同类物种只有 2 万这个样本的注释质量就值得怀疑后续结果解读也要慎重。4.3 物种文件命名与后续可读性OrthoFinder 本身对文件命名没有特殊要求但考虑后续 CAFE5 结果的可读性建议把文件名命名为“物种缩写 版本号”的形式例如Ath.fa Osa.fa Zmays.fa注意文件名不要包含点号以外的特殊字符更不要有()、%这类符号否则后续 shell 循环脚本里很容易踩坑。4.4 系统发育树文件准备CAFE5 运行需要一棵带分支长度的物种树ultrametric tree即所有叶子到根的距离相等。OrthoFinder 会输出一棵物种树SpeciesTree_rooted.txt但这棵树通常是“非超度量”的直接拿来跑 CAFE5 会报错或导致结果偏差。推荐做法是用r8s或PhyloBayes对 OrthoFinder 得到的物种树做时间校准得到超度量树。这里我不展开讲 r8s 的操作因为很多在线教程讲得比我详细。但如果你只是做初步分析也可以用 R 包ape的chronos()函数来实现时间校准代码极简效果也够用library(ape) tree - read.tree(SpeciesTree_rooted.txt) ultra_tree - chronos(tree) write.tree(ultra_tree, file SpeciesTree_ultrametric.nwk)这里输出的SpeciesTree_ultrametric.nwk就是后续 CAFE5 输入的树文件。5. 实操全程从 OrthoFinder 到 CAFE5 的完整命令流5.1 运行 OrthoFinder数据准备齐全后运行这条命令orthofinder -f proteomes/ -t 16 -a 16参数解释-f指定输入文件夹-t是 blast 用的线程数-a是分析流程用的线程数如果内存吃紧可以加上-M msa让多序列比对用 MSA 模式而不是默认的基因树模式但那会拖慢速度。跑完之后在proteomes/下会生成OrthoFinder/Results_日期/目录。我们关心的文件有三个cd proteomes/OrthoFinder/Results_*/Orthogroups/Orthogroups.tsv最核心的矩阵文件Orthogroups_UnassignedGenes.tsv没有被分到任何 orthogroup 的孤儿基因列表Orthogroups.txt每个 orthogroup 里的所有基因 ID用于后续提取序列。5.2 过滤单拷贝基因家族可选但推荐如果你后续还想做“基于单拷贝直系同源的物种树构建”或“共线性分析”可以在这一步把“所有物种都有且仅有一个拷贝”的基因家族挑出来。用下面这个 awk 命令awk -F\t NR1{for(i2;iNF;i) n[i]1; next} {ok1; for(i2;iNF;i){if($i!1) ok0} if(ok) print} Orthogroups.tsv Orthogroups_singlecopy.tsv即使不为了建树单拷贝基因家族也可以作为后续“看具体家族进化树”的候选集合建议保留一份。5.3 将 Orthogroups.tsv 转换为 CAFE5 输入CAFE5 输入文件需要包含三部分内容第一行物种数、基因家族数第二行物种 ID 列表用 Tab 分隔第三行开始家族 ID 每个物种的基因数量。我们可以直接用 Perl 脚本转换。我自己常用的脚本是orthogroups_to_cafe5.pl逻辑清楚且能处理文件头#!/usr/bin/perl use strict; use warnings; my $tsv shift ARGV; open(my $fh, , $tsv) or die Cannot open $tsv\n; my $header $fh; chomp $header; my sp split /\t/, $header; shift sp; # 去掉第一列 Orthogroup my $n_sp scalar sp; my rows; while (my $line $fh) { chomp $line; my cols split /\t/, $line; my $og shift cols; # 跳过全零的家族 my $sum 0; $sum $_ for cols; next if $sum 0; push rows, [$og, cols]; } my $n_fam scalar rows; print $n_sp\t$n_fam\n; print join(\t, sp), \n; for my $row (rows) { print join(\t, $row), \n; }保存脚本把Orthogroups.tsv喂进去就能得到 CAFE5 的输入文件perl orthogroups_to_cafe5.pl Orthogroups.tsv cafe5_input.txt排错提示如果cafe5_input.txt里的物种名和树文件里的物种名不一致CAFE5 会直接报错。所以树文件里的 tip 标签必须和这里的列名完全一致。这个坑我踩过不止一次务必核对。5.4 运行 CAFE5 并解析输出拿到输入文件和超度量树之后运行 CAFE5cafe5 -i cafe5_input.txt -t SpeciesTree_ultrametric.nwk -o cafe5_result/ -p 0.05参数说明-i输入文件-t树文件-o输出目录-p显著性阈值默认 0.05如果你的数据量很大还可以加-c指定线程数。运行结束后到输出目录看结果cd cafe5_result/重点看两个文件Base_family_results.txt每个基因家族的统计结果包括 family ID、家族-wide P-value、每个节点的家族大小估计等Base_change_results.txt每个节点上的扩张/收缩显著性结果。Base_change_results.txt长这样message: The p-value denotes the probability of the number of gains/losses... Node Family P-value Significant 12 OG0000000 0.001 Y 12 OG0000001 0.230 N ...这里的 Node 编号对应树文件中的内部节点编号。要搞清楚每个节点对应哪个分支我一般用FigTree打开带注释的树文件CAFE5 会生成Base_asr.tre等注释树或者用 R 包ggtree来可视化。如果你不想搞可视化至少要知道内部节点编号和分支的关系判断哪个物种谱系发生了显著变化。6. 数据后处理与结果可视化6.1 提取显著扩张/收缩的基因家族拿到Base_change_results.txt后筛选出显著变化的家族awk -F\t $4Y Base_change_results.txt | cut -f2 | sort -u significant_families.txt然后从原始Orthogroups.tsv里提取这些家族的详细成员信息方便后续富集分析。6.2 用 R 绘制家族大小变化热图我习惯用 R 的pheatmap对“物种 × 基因家族”矩阵画热图以直观展示哪些家族在特定物种里显著扩张。可以用Orthogroups.tsv做 log2 转化后画图library(pheatmap) data - read.delim(Orthogroups.tsv, row.names 1) # 过滤掉全零行 data - data[rowSums(data) 0, ] # 取变异最大的前 1000 个家族展示 data - data[order(apply(data, 1, var), decreasing TRUE)[1:1000], ] pheatmap(log2(data 1), show_rownames FALSE, scale row)这个热图非常适合放在文章补充材料里审稿人也爱看。6.3 GO/KEGG 富集分析对显著扩张的基因家族的成员基因做 GO 或 KEGG 富集分析是补充生物学故事的必要环节。具体工具取决于你的物种注释程度。如果物种有现成的 GO 注释比如从 Ensembl 下载用 R 包clusterProfiler就能跑如果注释不完善可以先用InterProScan做一个功能注释再走富集流程。这里不展开编码因为不同物种差异太大但有一个通用建议不要只拿显著扩张的家族做富集一定要有背景集。背景集一般用所有 orthogroup 覆盖的基因集否则结果会有严重偏差。7. 避坑指南这 6 个坑我踩过你别再踩了7.1 树文件不是超度量树CAFE5 最严格的要求就是树必须满足超度量性质ultrametric。如果你直接拿 OrthoFinder 输出的树跑 CAFE5很可能会遇到这类报错Initialization failed: tree is not ultrametric解决办法就是用前面提到的chronos或r8s做时间校准。顺便提醒一下如果你用chronos()它的默认平滑参数对大型树可能跑得很慢可以适当调整lambda值。7.2 输入矩阵里有“全零行”或极端异常值Orthogroups.tsv 里有一些家族可能在任何物种中都没有基因这听起来不可思议但确实会发生尤其是孤儿序列被独立聚成组后又被过滤的情况下。如果 CAFE5 遇到所有物种都为 0 的家族它会直接崩掉或者给出无意义的结果。所以转换脚本里的next if $sum 0那行过滤逻辑非常必要。7.3 细胞器基因和低复杂度序列没过滤线粒体、叶绿体基因以及大量低复杂度序列比如富含半胱氨酸的重复区域会干扰 BLAST 和聚类。如果你是做核基因组编码蛋白的分析建议先用seqkit或自定义脚本把细胞器来源的序列剔除至少也要检查你的输入 FASTA 里有没有线粒体/叶绿体来源的蛋白。7.4 物种命名不一致导致 CAFE5 报错这是最烦人但也是最常见的问题。OrthoFinder 运行时自动生成的物种树里物种名是“输入文件名去掉扩展名”得到的。而你在cafe5_input.txt表头里的物种名可能写的是全名或别名。两边一旦不一致CAFE5 会说树里有标签在数据中找不到。我的习惯是在运行 OrthoFinder 之前就把文件名定好所有后续步骤都用同一套物种名绝不中途改名。7.5 只盯着 P-value忽略了分支长度的影响CAFE5 的推断高度依赖树的分支长度。如果一棵树的分支长度校准得奇差无比比如某条分支比所有其他分支长 10 倍那 CAFE5 大概率会把大多数家族都判定为显著变化这种结果是不可信的。所以我强烈建议在跑 CAFE5 之前用FigTree打开你的超度量树看一眼如果存在明显异常的长分支优先回到物种树构建步骤排查问题而不是盲目继续。7.6 忽略了物种树的拓扑不确定性基因家族收缩扩张分析的结果会被物种树拓扑影响。如果你的物种树在某些分支上支持率很低比如只有 40% 的 bootstrap 支持率那这些节点上的扩张/收缩推断结果就缺乏可靠性。吧里有些人的做法是“跑多棵物种树然后取一致结果”这个思路成本较高但严谨度确实更高。至少你应该在文章里说明树是用什么方法构建的、支持率如何。8. 实操心得与后续扩展方向整套流程跑下来我最大的感受是OrthoFinder CAFE5 的组合看起来只有两个软件真正耗时间的不是运行而是前面的数据清洗和后面的结果解读。OrthoFinder 本身跑得非常快5 个物种的蛋白组16 线程大概 1~2 小时就结束了。CAFE5 如果做全基因组范围的分析也就是几分钟到十几分钟的事情。相比之下我花在过滤序列、校准树、核对物种名上的时间至少是运行时间的 3 倍。另外多说一句关于“显著扩张/收缩家族”的生物学解读。我一直提醒学生CAFE5 是基于序列计数统计出的“显著变化”它不等同于“该基因家族的扩张就是这个物种适应性演化的原因”。要支撑这样的结论你还得补充表达量数据、选择压力分析dN/dS、结构域保守性分析甚至功能实验。但在基因组学的初步探索阶段CAFE5 的结果已经足以提供非常有价值的候选基因列表再结合人工筛查完全可以支撑起一篇高质量的比较基因组学文章。扩展方向上如果你对某个感兴趣的具体基因家族比如 NLR 抗病基因家族、细胞色素 P450 家族做细致分析可以在 OrthoFinder 结果里抽取该家族的序列用 MAFFT 比对后用 IQ-TREE 构建基因树再结合物种树做基因树-物种树一致性分析比如用 NOTUNG 或 R 的 phytools 做 reconcile。这套组合拳基本覆盖了基因家族进化的主流分析框架。最后分享一个脚本层面上的小技巧在跑任何一步大型分析时都顺手记录一下软件版本、参数和运行时间最好写进一个analysis_log.txt里。因为后面写作时你一定会被要求提供“详细分析方法”部分到时候现回忆参数是用 0.05 还是 0.01、树是用 chronos 还是 r8s准得头疼。我是吃过这个亏的现在每次跑流程都开一个 log写文章时直接复制粘贴就能用省事太多。以上就是我从 OrthoFinder 跑到 CAFE5 的完整流程和踩坑经验。照着做你大概率能少走一半弯路。
