企业官网建设流程全解析

做比较基因组学项目，基因家族鉴定和注释跑完之后，最常被审稿人追问的一件事就是：你这些家族在物种间的拷贝数差异，到底对应什么生物学事件。基因家族扩张收缩分析（CAFE5）就是回答这个问题的标准动作。这篇文章我把从输入数据准备、参数调整、output解析到最终可视化的完整流程写出来，包含我在实际项目中踩过的坑和总结的绘图习惯。适合正在做物种基因组进化的同学参考，也适合刚拿到OrthoFinder结果、打算用CAFE5出图的初学者直接抄作业。

1. 为什么CAFE5是"家族变化"这个问题绕不开的工具

1.1 扩张收缩不是锦上添花，而是物种适应故事的核心证据

很多人会有个误区，觉得基因家族分析做到注释和GO/KEGG富集就够了。实际上，如果一篇比较基因组学文章里没有"哪些基因家族发生了谱系特异性扩张、哪些收缩"这部分内容，审稿人几乎一定会问：家族层面的变化是怎么样的？原因很简单，基因家族的拷贝数变化是表型进化的一个直接线索——比如植食性昆虫里解毒酶家族的扩张、深海环境里视蛋白家族的收缩，都是这一类分析给出的强证据。

CAFE全称是Computational Analysis of gene Family Evolution，第五代版本（CAFE5）是目前应用最广的实现。它的核心逻辑可以理解成：你有一棵物种树，还有每个物种里每个基因家族的基因数，CAFE5用一个随机的出生-死亡过程去模拟家族大小在树枝上的演化，估计出每个家族在每条谱系中的获得/丢失速率，然后判断哪些分支上的家族大小变化显著偏离随机期望。用一句大白话讲，它就是在问：给我一棵树和一些家族大小数据，哪些地方的变化"大得不像是碰巧发生的"。

1.2 CAFE5相比老版本到底改进了什么

如果你在服务器上还留着CAFE 2.x或3.x的老脚本，我建议尽快换。CAFE5从4.x开始最大的变化就是支持多线程并行，我以前跑几十个物种、几万个家族的矩阵，老版本可能要跑一晚上，现在CAFE5指定-c 8之后，半小时到一小时基本能出结果。此外，CAFE5对输入矩阵的容错更友好，还能直接估计多lambda模型——也就是允许不同分支拥有不同的演化速率，而不是全树只用一个速率。

另一个值得提的点是错误率模型（error model）。基因家族计数矩阵是从OrthoFinder之类的软件得到的，这个计数本身是有误差的，比如基因组组装质量差会导致某些家族计数偏低。CAFE5提供的错误率模型会在似然计算中把这类误差考虑进去，代价是计算量明显上涨。我的建议是：如果项目只是初步看趋势，先不开error model，把候选家族筛出来；如果到写文章阶段，再用带error model的结果做最终确认。

2. 跑通CAFE5的最低配置：输入文件长什么样、参数怎么定

2.1 基因家族计数矩阵的格式与生成方式

先看最关键的文件：计数矩阵。CAFE5的输入并不是常见的行列矩阵格式，而是三列的记录形式。

Family_id Species_name Gene_count Family_001 Oryza_sativa 5 Family_001 Zea_mays 4 Family_001 Arabidopsis_thaliana 2 Family_001 Glycine_max 8 Family_001 Solanum_lycopersicum 3 Family_001 Sorghum_bicolor 6

如果你用的是OrthoFinder结果，比较省事的方式是用它自带的转换脚本，把Orthogroups.tsv转成CAFE5的计数文件。

python OGC2CAFE.py -i Orthogroups.tsv -g species_abbreviation.txt -o cafe_input.txt

这个脚本生成的cafe_input.txt就已经是标准格式。如果你是自己写脚本整理，要格外注意：每个family在不同物种中都要出现，没有基因就写0，不要留空；物种名列不要有多余空格或_前后不一致的情况。我遇到过的最多的报错就是物种名这里对不上，导致后面解析树文件时直接退出。

2.2 带枝长的超度量树：最容易忽略的前提

树文件的要求比很多人想的严格。CAFE5要求输入一个带枝长的Newick树，而且最好是超度量树（ultrametric tree），也就是从根到所有叶子的遗传距离相等。这个要求背后的原因是：CAFE5模型假设所有现有物种经历了相同的演化时间总量，枝长只是用来把进化事件分配到不同段落上。如果你的树没有经过时间校准，根到不同物种的枝长总和差很多，算出来的速率估计就会失真。

检查是否超度量的方法很简单：

library(ape) tree <- read.tree("species_tree.nwk") is.ultrametric(tree)

如果返回FALSE，不要直接丢给CAFE5跑，先做时间校准。有化石校准信息的话用treePL或r8s，没有的话至少用chronos这种基于分子数据的方法做个相对时间树。如果只是快速验证流程，也可以用phytools::force.ultrametric强制转化，但文章里不建议这么干。树文件格式大概是这样：

(Oryza_sativa:0.12,(Zea_mays:0.08,(Sorghum_bicolor:0.06,(Glycine_max:0.10,(Arabidopsis_thaliana:0.05,Solanum_lycopersicum:0.05):0.08):0.04):0.05):0.03);

2.3 参数取舍：-p、-k、-l、-c分别意味着什么

CAFE5的命令行参数并不多，但每个都影响结果。以我常用的命令为例：

cafe5 -i cafe_input.txt -t species_tree.nwk -o cafe5_out -p 0.01 -k 2 -l -c 8

-p是初始lambda值，也就是出生-死亡速率估计的起始猜测。这个值不必追求精确，一般取0.001到0.1之间都行。有一点要注意：CAFE5的结果对起始值有一定敏感性，所以严谨的做法是换几个不同的-p值跑一遍，看看显著家族列表是否基本一致。如果一个家族在某个种子下P值0.01、另一个种子下P值0.3，那它基本处于边缘状态，写文章时我会选择不放图里。

-k是最重要的模型选择参数，表示允许的lambda类别数。-k 1是全树用一个速率，#最稳妥但也最不灵活；-k 2以上是假设不同分支可以有不同的速率，计算量会涨，但更贴近真实演化场景。我个人的经验是先跑一次-k 2 -l看看整体似然是否明显提升，再用似然比检验决定是否采用多速率模型。-l控制是否用最大似然法搜索最优lambda，建议常规分析都加上。-c是线程数，服务器内存够就尽量给大，比如8或16，能显著缩短全家族扫描时间。

3. 从输出文件到可视化数据：中间格式到底该怎么读

3.1 Base_results.txt：先看全貌再谈绘图

CAFE5的输出目录里第一个要打开的文件是Base_results.txt（实际文件名前缀取决于是否用了--cores、--lambda等选项，但结构一致）。它是一个扁平的大表格，每行对应一个基因家族，列包括家族号、P值、估计lambda、是否显著等。

Family_id pvalue lambda significance Family_001 0.001 0.003 Yes Family_002 0.42 0.001 No Family_003 0.00005 0.001 Yes

拿到这个文件后，第一件事不是急着画图，而是先统计显著家族的数量。我一般会在R里快速过一遍：

df <- read.delim("Base_results.txt") summary(df$pvalue) sum(df$pvalue < 0.05)

如果显著家族数量特别多（比如超过家族总数的30%），通常说明模型设定有问题，比如树不满足超度量、计数矩阵里的0过多导致分布极端，或者你更需要一个带error model的结果。如果显著家族太少甚至没有，则要考虑是不是过滤阈值太严，把很多低频家族删掉了。

3.2 Base_change.tre：把变化事件直接映射到系统发育树上

这个文件是整个可视化流程的核心，它本质上是一棵Newick树，但每个节点的标签里嵌入了家族大小变化信息。例如某个节点可能显示为Oryza_sativa[-4]或[+3]Node12，具体格式会随CAFE5版本略有变化，但思路一致：方括号里的数字表示该分支上总共发生了多少次扩张（+）或收缩（-）事件。

读取它的方式和读普通树完全一样，关键是在绘图时把标签解析成可用的数值。我通常用一条正则把方括号里的内容抽出来：

library(ggtree) library(phytools) change_tree <- read.tree("Base_change.tre") # 提取尖括号内的变化标记 get_change <- function(label) { m <- regmatches(label, gregexpr("\\[([+-]?[0-9]+)\\]", label))[[1]] if(length(m) == 0) return(0) as.numeric(gsub("\\[|\\]", "", m[1])) }

这一步做完之后，这棵树的每个节点都对应一个"变化强度"的数值，接下来所有统计图和树图都建立在它上面。

3.3 Base_families.tsv：从全家族结果中筛选真正有故事的基因

Base_families.tsv记录了每个家族在每个节点上的具体推断值，是后续深入分析（比如看某个扩张家族在特定物种里的拷贝数变化轨迹）的基础。不过做可视化之前，更常用的是它包含的P值和家族ID，用来和Base_results.txt对应。我的筛选流程是这样的：先按P值取显著家族，再做BH多重检验校正，最后只看FDR < 0.05的。

df <- read.delim("Base_families.tsv") df$fdr <- p.adjust(df$pvalue, method = "BH") sig_families <- df[df$fdr < 0.05 & df$pvalue < 0.05, ]

筛选后建议把显著家族分成扩张和收缩两类。判断标准可以看家族在现有物种里的平均拷贝数相比祖先状态是增加还是减少——CAFE5输出的祖先状态就存在Base_asr.tre里，也可以用更粗略的方式：某个家族在这个谱系的末端物种中基因数中位数 > 所有物种中位数的一半以上，就归为扩张方向。文章里通常只挑选排名前几的代表性家族来画图，因为所有显著家族全画出来，图会非常挤。

4. 四张必出图详解：从树到热图，每一张图怎么画

4.1 树图：在分支上标注扩张和收缩事件数

这是基因家族扩张收缩分析最经典的一张图，几乎每篇比较基因组学文章里都会出现。它把Base_change.tre读进来，用ggtree画系统发育树，然后用节点标签颜色或者点大小表示扩张/收缩强度。画法如下：

library(ggtree) library(ggplot2) p <- ggtree(change_tree, size = 0.8) + geom_tiplab(size = 3, offset = 0.01) + geom_nodelab(aes(label = label, color = as.numeric(label)), hjust = -0.3, size = 3) + scale_color_gradient2(low = "blue", mid = "white", high = "red") + theme(legend.position = "right") ggsave("family_change_tree.pdf", p, width = 8, height = 6)

这里geom_nodelab直接使用节点标签数值做颜色映射。如果是第一次画，我建议先只标数字、不上色，确认解析结果没问题再做颜色映射，因为Base_change.tre的标签格式在不同版本之间存在差异，盲目画很容易出现所有节点都是同一个颜色。另外，图里的分支粗细建议用"事件数量"映射，因为大多数读者看到粗分支的第一反应是"这里变化大"。

4.2 统计柱状图：不同物种分支的扩张和收缩数量对比

树图适合看全局，但如果要强调某个分支比其他分支发生了更多扩张或收缩，柱状图更直观。做法是从Base_change.tre中提取每个末端物种的变化量，然后按扩张/收缩分组：

library(dplyr) library(ggplot2) # change_data 包含 species, type, count 三列 p <- ggplot(change_data, aes(reorder(species, count), count, fill = type)) + geom_col(position = position_dodge(width = 0.7), width = 0.6) + coord_flip() + labs(x = NULL, y = "Gene family count") + theme_minimal()

实测下来，柱状图配合树图放在同一个figure里效果最好：上边是树，下面是对应的柱状图，两者共享物种顺序，读者一眼就能把"哪个分支变化大"和"具体变化了哪些家族"对应起来。这种情况适合投文章主图，也适合做组会的汇报图。

4.3 饼图或环形图：整体扩张收缩比例

这张图适合放在补充材料里，用来快速交代全基因组的家族变化格局。统计逻辑很简单：在Base_change.tre中，把每个节点标签里的正负号分别累加，得到总扩张事件数和总收缩事件数，然后画饼图。严格意义上讲，"家族数量"和"事件数量"是两回事——一个显著扩张家族在某个分支上可能贡献了多次事件，但很多文章直接用事件数量代表家族数量，我建议你在图注里写清楚统计口径。

饼图我一般会在R里用一行代码：

pie(table(event_type), col = c("#E64B35", "#4DBBD5"))

如果觉得base R的饼图太朴素，可以转成ggplot2的柱状图加coord_polar()。但无论怎么画，这张图信息量有限，放在正文里有点浪费版面，除非你们的审稿人明确要求看全基因组统计比例。

4.4 热图：显著家族在多个物种中的拷贝数模式

最后一张是针对具体显著家族的细节图。做法是取筛选出的显著家族，把这些家族在所有物种中的原始基因数整理成矩阵，做个带聚类和注释的热图。我的做法是行的顺序按系统发育树排序，列按家族聚类，颜色用log2(基因数+1)转换，避免个别家族拷贝数特别大拉高色阶。

library(pheatmap) # sig_matrix: 行是家族，列是物种，值是基因数 pheatmap(log2(sig_matrix + 1), scale = "row", cluster_rows = TRUE, cluster_cols = FALSE, show_rownames = FALSE, treeheight_row = 20)

这张图的价值在于展示"同一个家族在不同物种中拷贝数的差异"。比如某家族在水稻里有20个拷贝，在拟南芥里只有2个，对应到树图上就是水稻这一支的强扩张信号。如果你有条件，还可以在热图旁边标注每个家族的注释信息，比如激酶、抗病基因、转录因子这些功能类别，会提高整张图的信息密度。

5. 实操中容易翻车的四个细节：排查与修正

5.1 树文件不满足超度量导致报错或结果失真

报错一般是这种情况：CAFE5运行到中途提示树的枝长总和不等，或者干脆在读取树文件后立即报ERROR: tree has different root-to-tip distances。这个问题我在早期项目里碰到过很多次，因为直接从OrthoFinder或已有文献里拿到的树，并不会自动满足超度量要求。

排查路径是这样的：先在R里用is.ultrametric()检查，如果是FALSE，再算一下每片叶子的root-to-tip距离到底差多少。如果差异很小（比如千分之一），你可以用phytools::force.ultrametric()强制转换；如果差异很大，说明树本身不是一棵时间树，需要回到分子序列或化石校准流程重做时间树。实在只能手动调枝长的话，记得在图例里注明原始树来源。

5.2 计数矩阵零值过多导致计算量暴涨

CAFE5在做全家族扫描时，完全全部为0或全部为1的家族，其对似然的贡献是恒定的，但计算过程仍会消耗大量时间。比较稳妥的做法是在跑CAFE5之前做一次过滤：

• 在所有物种中基因数均为0的家族直接删除；
• 单拷贝在所有物种（都是1）的家族，理论上不会产生任何扩张收缩信号，通常也建议删掉；
• 设定一个最低非零物种数阈值，比如至少5%的物种里基因数大于1，否则过滤掉。

这个步骤能让你少跑近一半的无效家族，对大型基因组项目来说省下的时间非常可观。需要注意是：如果用--error_model，不要过滤太狠，因为错误率模型需要足够的数据量来估计计数误差。

5.3 P值全显著或全不显著的"两头极端"问题

如果跑完结果发现所有家族P值都小于0.05，先别高兴太早，这基本说明模型有问题，而不是你运气好。最常见的原因是零膨胀——计数矩阵里有大量0，使所有家族都被推向了极端分布，此时应该有条件地采用误差模型。另一个原因是树枝长没有正确的绝对时间尺度。反过来，全部不显著时，先检查是不是过滤时把变异大的家族也删掉了。

我一般的处理顺序是：先重新看计数矩阵的分布直方图（绝大多数家族是0还是1个拷贝），再把过滤阈值放宽，最后增加lambda类别数（从-k 1换成-k 2）。多数情况下，"全是显著"或"全不显著"都回落到一个合理的中间状态——大约5%-15%的家族显著。

5.4 结果不稳定：随机种子与多次运行的取舍

CAFE5的-p参数内部包含随机数初始化，所以严格意义上结果会受启动种子影响。我在实际项目里会挑三个不同的初始 lambda（比如 0.001、0.01、0.1）各跑一遍，然后把三份Base_families.tsv合并起来比较。如果某家族在两次以上运行中都显著，基本可以放心收入正文；如果只有一次显著，多半是处于边缘状态。这个操作看似费时间，但能避免文章修稿时因为复现结果不稳定被审稿人挑战。

另外，如果你有强背景知识认为某个分支演化速率极高或极低，可以显式地通过-k 2指定两套速率类别，再在结果里看该分支被分配到哪个速率。不要把所有的判断都丢给软件，领域先验知识在这个环节很重要。

6. 写到最后：关于基团家族结果解释的一点体会

CAFE5跑完、图也画完了，并不意味着分析就结束了。我这两年最深的感受是：一张漂亮的扩张收缩图只是第一步，真正的问题在于"为什么"。如果一个抗病相关基因家族在某个谱系里显著扩张，你是否能回到基因组组装和注释里确认这些拷贝不是碎片或者伪基因？如果一个转录因子家族显著收缩，你是否能从表达数据里看到对应靶基因的补偿性变化？这些后续验证往往决定了这篇文章能不能从"组学描述"上升到"生物学发现"的层次。

在工具链上再啰嗦一句：CAFE5输出的是点估计和P值，但基因家族演化本质是随机过程，对单个家族的结果解释要克制。我会在文章里给每个显著家族加一行"这个家族是否在多次独立分析、不同方法（比如OrthoFinder和InParanoid）中表现一致"的内容，这也是外审比较看重的一个严谨性信号。分析本身不难，难的是你想要讲一个什么样的演化故事，CAFE5和它的可视化只是帮你把故事的第一页写出来而已。