CAFE5基因家族扩张收缩分析：从原理到实操全指南

1. 基因家族扩张收缩分析，到底在解决什么问题

做基因组项目的人，迟早会撞上一个问题：手上有了一个刚组装、注释完的基因组，接下来除了做共线性、做进化树，还能做什么？如果这个物种恰好有近缘物种的基因组可以参考，那基因家族扩张收缩分析基本是绕不开的一步。

这个分析的核心逻辑并不复杂。同一个基因家族在不同物种里，成员数量往往不一样。有些家族在某个物种里特别“膨胀”，比如与抗病相关的NBS-LRR类基因在不少植物里动辄几百个成员；有些家族则“萎缩”得厉害，甚至彻底丢失。这种数量上的差异背后，往往藏着适应性演化的线索——某个家族显著扩张，可能是因为环境压力、共生关系或者特殊的生理需求在驱动。

CAFE5就是专门干这件事的工具。它全称是Computational Analysis of gene Family Evolution，目前最新版本是第5代。它做的事情，是在一个给定的物种系统发育树上，根据每个基因家族当前的成员数量，反向推断祖先节点的家族大小，然后通过统计检验找出哪些分支上的哪些家族发生了显著扩张或收缩。

这里有个容易混淆的点，必须说清楚：CAFE5分析的“扩张/收缩”不是简单比较两个物种的家族成员数谁多谁少，而是沿着系统发育树，在进化时间尺度上建模，找出在特定分支上偏离了背景速率的家族。换句话说，它比“你多我少”要严谨得多。

这个分析适合谁来用？只要你的课题里有一个刚测序的基因组，并且有两三个近缘物种的蛋白序列可以做基因家族聚类，就可以跑。不需要特别深的数学背景，但需要一点耐心去准备输入文件。我个人觉得，CAFE5是基因组进化学分析里性价比很高的一步——它花不了多少计算资源，却能在文章里贡献一张很有说服力的图。

2. 跑CAFE5之前的准备工作，卡住大多数人的地方在这儿

CAFE5本身跑起来很快，真正让人头疼的是输入数据的准备。这一步做好，后面基本就是流水线操作。

2.1 输入文件之一：带分支长度的系统发育树

CAFE5需要的树不是那种只有拓扑结构的树，而是一棵带有分支长度（通常代表时间或替换数）的树，单位是时间或者替换数，一般用基因树的替换数来近似。

我当时用的是OrthoFinder跑出来的物种树，文件格式是newick，类似这样：

(((speciesA:0.12,speciesB:0.15):0.08,speciesC:0.20):0.1,speciesD:0.24);

有几个坑要先提醒:

• CAFE5对树的分支长度单位没有硬性要求，但分支长度如果过小，比如低于0.01，可能会出现计算上的数值问题。我自己遇到过几次“Error: branch length too small”的情况，后来统一把分支长度乘以100，问题就解决了。本质上缩放不影响结果，因为似然计算的分支长度是相对的。
• 树的物种名必须和后面计数矩阵的列名严格一致，连大小写、下划线都不能差，否则会直接报“Species not found”。
• 如果树没有根，CAFE5会默认在最长分支的中点找根，但我建议自己先给定好根的OuGroup。

2.2 输入文件之二：基因家族计数矩阵

这个矩阵就是每个物种里每个基因家族的成员数量。行是基因家族，列是物种，格式很直观：

FamilyID speciesA speciesB speciesC speciesD OG0000001 12 8 15 5 OG0000002 3 4 2 6

这个矩阵怎么来？我通常的做法有两种：

第一种是用OrthoFinder的结果做清洗。OrthoFinder会输出一个Orthogroups.tsv文件，每一行是一个直系同源群，列是物种，单元格里是该物种在这个家族里的所有基因ID。我写了一个小脚本，把每列按逗号分隔后统计个数，直接生成计数矩阵。这里要注意，OrthoFinder里的孤儿基因（即单拷贝至多拷贝的独有基因）如何处理需要自己拿捏。如果某个家族只有一个物种有基因，其他物种都是0，这种家族建议过滤掉，因为它在似然计算中提供不了什么有效信号，反而可能引起歧义。

第二种是用InterProScan的结果自己做家族定义。比如你对某个特定结构域家族感兴趣，可以把你关注的基因ID列表整理成家族，再统计每个物种的基因数。这种做法的好处是分析更聚焦，缺点是需要额外做一次domain注释，工作量略大。

2.3 过滤标准：别把所有家族都扔进去

CAFE5官方推荐的做法是过滤掉成员数过多或过少的家族。我个人的经验如下：

• 成员总数为0的家族必须去掉。
• 所有物种里成员总数之和超过100的家族，建议剔除。因为家族太大，计算量会成倍增加，而且这些超大家族本身的演化模式就很特殊，容易干扰整体分析。
• 某一行里，如果有一个物种缺少该家族（计数为0），这个没关系。但如果该家族在超过一半的物种里都是0，我倾向于直接过滤，以减少不确定性。

实际跑下来，过滤后大概剩下一万多个家族是比较常见的情况。过滤条件可以写进一个小脚本，按行判断，几行代码就搞定。

准备完这两个输入文件，CAFE5的“原料”就够了。

3. CAFE5安装与运行实操，附命令行模板

3.1 安装：两种方式，建议直接conda

CAFE5的安装不算复杂，两种主流方式：

用conda的话一句话搞定：

conda install -c bioconda cafe5

如果conda源有问题，也可以从GitHub源码编译。依赖主要是libsqlite3、libmysqlclient这些，编译过程也比较顺利，但需要手动指定一些库路径，没有conda省心。我自己在服务器上装的时候遇到过一次编译报错，原因是libsqlite3版本过低，后来用conda重新建了个环境，几分钟就解决了。所以我的建议是：直接用conda建一个干净环境，别在主环境里折腾。

3.2 运行命令：基础用法和参数调优

CAFE5的命令行参数不算多，核心就这几个，一个最基础的运行命令长这样：

cafe5 \ -i ../input/family_counts.txt \ -t ../input/species_tree.txt \ -o ../output/cafe5_run1 \ -p 0.05 \ -c 8

参数的含义逐一说下：

• -i：输入计数矩阵文件。
• -t：输入物种树文件，newick格式，带分支长度。
• -o：输出目录，CAFE5会为每次运行创建一个新目录，不建议重复使用同一个目录。
• -p：p值阈值，用于筛选显著扩张/收缩的家族，默认0.05，也可以设成0.01更严格。
• -c：CPU核心数。CAFE5支持多线程，但实测下来线程数设太高时提升有限，取8-16即可。
• -k：这个参数很关键，它指定λ值的个数。默认是2，即让程序自动选择两个λ值来拟合数据。λ是基因家族的获得/丢失速率，CAFE会同时估算几个λ再比较哪个模型更合适。如果你不想花时间调优，就用默认值；但如果你观察到结果中的λ值分布不理想，可以尝试-k 1，让整个树共享一个λ，有时反而更稳定。

CAFE5运行速度很快，一万多个家族，8个线程，基本十几分钟就能跑完。当我第一次跑完看到输出目录里各种文件时，还是感慨这工具效率确实高。

3.3 一个要注意的细节：随机种子

CAFE5里面有一个计算步骤涉及随机性（在估算祖先状态时使用了随机化），所以如果你想让结果可重复，需要设置固定的随机种子，参数是-r。

比如-r 12345，如果不设，两次运行的数值会有微小差异，对主要结论影响不大，但审稿人如果要你提供精确复现流程，还是设一下比较好。

4. 结果文件逐个解读，CAFE5输出里到底有什么

CAFE5运行完之后，会在输出目录里生成一系列文件，第一次接触的人容易看得一脸懵。这里我把关键文件整理一下，按重要程度排序。

4.1 核心文件速查表

其中最常用的就是Base_change.tab和Base_family_results.txt。

Base_family_results.txt每一行是一个家族，最后几列会有p-value和db（方向和显著性标记），你可以用awk命令直接筛出显著家族：

awk '$NF < 0.05' Base_family_results.txt | wc -l

Base_change.tab的列结构是：家族ID + 每个节点的家族大小变化，正值代表扩张，负值代表收缩，用这个文件可以精确追溯一个家族在每个内节点上的演化过程。

4.2 快速定位显著扩张/收缩家族

我一般拿到结果后会做一个这样的筛选：

1. 先从Base_family_results.txt里筛出p<0.05的家族；
2. 再从Base_change.tab中找到目标分支（比如你关注的某个谱系）上有显著变化的家族；
3. 将这些家族成员基因ID提取出来，做GO/KEGG富集分析。

这个流程走下来基本能锁定几个候选家族，后期再用蛋白结构或表达数据进一步验证，就是很典型的套路了。

5. 可视化这一步，把进化故事讲清楚

CAFE5跑完只是数据分析的中间站，怎么把结果转化成论文里直观的图，才算真正完成闭环。

5.1 用R语言绘制扩张收缩数量分布图

最常见的可视化方式，是在系统发育树旁边画每个分支上扩张和收缩家族的数量条形图。这一步我用的工具是ggtree+ggplot2。

先提取每个分支上的显著扩张/收缩家族数量，整理成这样的格式：

Branch Expansion Contraction A 156 34 B 89 101 C 12 56

然后用ggtree读取物种树，加上geom_text和geom_bar，调整好色板之后出图。这种图的特点是信息量集中在一张图里，能一眼看出哪个谱系扩张最猛、哪个谱系收缩最多，非常适合放在文章主图里。

5.2 演化速率可视化：换个角度挖掘信息

除了数量和p值，CAFE5还会输出每个家族的λ值（获得/丢失速率），这个信息往往被忽略。我最近一次分析中尝试把每个家族的λ值提取出来，和家族大小做散点图，发现一个规律：家族越大，λ越小，两者呈明显的负相关。这个图非常适合放在文章补充材料里，作为模型合理性的佐证。

还有一种可视化角度：热图。把显著扩张/收缩家族在不同物种间的成员数做成热图，结合物种树聚类，能快速识别出某些家族在某个谱系中的一致性扩张模式。做法也不复杂，筛选出显著变化的家族后，提取原始计数矩阵的子集，用pheatmap画热图即可。

5.3 如果有R基础，试试这个脚本思路

这里给一个简化的R代码思路，方便大家按自己的数据结构调整：

library(ggtree) library(ggplot2) # 读取树文件 tree <- read.tree("species_tree.nwk") # 读取分支变化统计 change_data <- read.table("branch_changes.txt", header = TRUE) # 用ggtree画树，在尖端加上条形图 p <- ggtree(tree) %<+% change_data + geom_tiplab() + geom_bar(aes(x = Expansion, fill = "Expansion"), stat = "identity", width = 0.5) + scale_fill_manual(values = c("Expansion" = "#E64B35", "Contraction" = "#4DBBD5"))

需要注意的是，树尖端的物种顺序要和条形图的行顺序对应，否则图会乱。%<+%是ggtree提供的一种关联数据的方式，它会自动按tip label匹配，一般不用操心顺序问题。

5.4 工具选型：不是越复杂越好

做可视化的时候，有人喜欢直接上iTOL，有人喜欢全套R代码。我的观点是：如果只是画一个最终的展示图，iTOL确实方便，拖拽上传就能出图，而且在线的交互功能很友好。但如果要批量处理多种方案、多次修改，R脚本的可复现性高得多。我在实际项目中通常先用iTOL快速预览，确定展示样式后再用R脚本定稿出图。

另外提一句，CAFE5官网提供了一个Python脚本cafe5_draw_tree.py，可以快速画一个带家族大小分布的基础树图，给那些不熟悉R的人一个兜底方案。不过这个脚本出的图样式比较简单，适合初筛，不太适合直接作为发表级图片。

6. 我在实际项目中踩过的五个坑

这部分是我最想写的，因为CAFE5本身不难跑，真正耽误时间的往往是这些边边角角的问题。

6.1 坑一：分支长度过小导致报错

有一次我在处理一个分支长度在0.001级别的树时，CAFE5直接报错退出，报错信息大致是“branch length too small”。网上搜了一圈，最后我干脆把整棵树的长度等比放大100倍，问题迎刃而解。原因是CAFE5内部计算时可能对极小值做了数值下溢保护，等比缩放不影响相对长度，也就不会改变推断结果。

6.2 坑二：计数矩阵物种顺序和树不一致

CAFE5在解析时要求计数矩阵的列名与树的tip label对应，但顺序没有硬性规定。话虽如此，有一次我矩阵里的物种名带了拼接的版本号，比如“A_genes”，而树里是“A”，CAFE5直接全报错。检查了半天才发现是名字不匹配。所以预处理时一定要统一好命名规则，物种名保持纯ID，不要夹带注释信息。

6.3 坑三：超大家族拖慢分析

起初我把所有家族都扔进去，结果运行时间一下子从十几分钟变成了两三个小时。后来按官方建议过滤掉总成员数超过100的家族，速度立刻恢复正常。建议大家一开始就加上过滤步骤，别等到跑出事来再回头改。

6.4 坑四：祖先状态重建结果的解读陷阱

CAFE5输出的祖先状态只是基于当前树和家族大小的推断，并不是实验验证的结果，更不代表真实的祖先一定具有那么多的基因家族成员。我看过一些新手直接把祖先状态当成“事实”来写，这是有风险的。合理的说法应该是“CAFE5推断该节点祖先可能含有X个家族成员”，语气要严谨。

6.5 坑五：p值趋近于0的结果要谨慎对待

筛选显著家族时，很多家族p值直接是0（实际上极小）。这种超显著结果有两种可能：一是该家族确实在某个分支上经历了爆发式扩增，二是数据本身存在偏差，比如注释质量问题导致某个物种的基因数虚高。处理方法是对这些极显著家族做一个手动核查，挑几个基因做PCR或看转录组表达量佐证。

6.6 一个额外的小建议：保存完整运行日志

CAFE5运行过程中，日志文件Cafe5_log.txt记录了命令行参数、运行时间和λ估计值。分析做完之后，把这日志和输入文件、输出文件归档在一起，形成一套完整可复现的目录结构，对后续修改或写文章方法部分都会省很多事。

7. 常见问题速查，建议直接收藏

8. 怎么把这套分析用在你的文章里

如果你手头有一个新组装的基因组，CAFE5的扩张收缩分析可以这样设计进故事线：

第一步，和2-3个近缘物种做OrthoFinder聚类，得到基因家族；第二步，运行CAFE5，识别出目标谱系上显著扩张和收缩的家族；第三步，对这些家族做功能富集分析，找到富集的GO条目或KEGG通路；第四步，结合转录组或表型数据，验证这些家族是否真的在特定组织或条件下有异常表达。

在我做过的一个植物基因组项目中，发现某个物种的抗病相关家族显著扩张，而且这些基因在根组织中高表达。结合该物种生长在病害高发的环境背景，这个结果就很自然地成了故事的核心亮点。审稿人对这种“数据-功能-表型”的闭环通常评价不错。

最后再分享一个小技巧：CAFE5的结果不要只盯着一张总图，试着把显著家族单独拎出来和近缘物种做一次多序列比对，看看扩张出来的基因是否保留了关键功能结构域。有时候你会发现，所谓扩张其实是假基因化后的“数量增长”，这种细节会让你的分析质量上一个台阶。

实际跑过几次CAFE5之后我的体会是，这套分析本身不难，真正的功夫在前期的数据清洗和后期的结果解读。工具再强也只是辅助，能把进化故事讲合理、讲扎实，才是做这个分析最值钱的地方。