从OrthoFinder到iTOL：基因家族进化分析全流程实战指南

做比较基因组学或者基因家族进化分析的人十有八九都卡过同一道坎手里的基因组和转录组数据堆了一堆想理清楚某个基因家族的演化历史却不知道从哪一步下手。OrthoFinder能帮你把多物种的蛋白序列聚成直系同源组iTOL能把基因树画得发表级好看但网上讲OrthoFinder的教程多半只讲到输出一堆文件夹就没了讲iTOL的教程又默认你手里已经有树文件。真正让人头疼的“中间环节”——怎么从OrthoFinder结果里精准捞出目标基因树、怎么处理序列比对和建树参数、怎么让iTOL的注释文件一次上传就通过——反而没人系统讲。这篇文章我就以兰科植物为例把这套“OrthoFinder → 目标OG筛选 → 多序列比对 → IQ-TREE建树 → iTOL美化”的全流程完整走一遍。兰科是被子植物里物种多样性最夸张的科之一蝴蝶兰、石斛、天麻这些大家都很熟做花发育、物种适应相关基因家族的演化分析时这套流程几乎必用。文章里所有命令都直接给了可复制的版本参数我也逐个解释为什么这么设。无论你是刚入门生信的学生还是已经跑过几次OrthoFinder但没摸透结果文件的老手照着做都能少踩几个坑。1. 流程整体设计为什么是OrthoFinder到iTOL这条链路1.1 这套流程能回答什么科学问题你需要先想清楚做基因树到底是为了回答什么问题最常见的两类第一是基因家族演化分析比如某个转录因子家族在兰科里发生了哪些扩张和收缩某个谱系特有分支是怎么来的第二是功能验证前的进化背景评估比如你克隆了一个新基因想知道它属于哪个亚家族、和已知功能基因的亲缘关系如何。这两种需求都跑不出同一个底层的分析链路先把多物种的所有蛋白序列按同源关系分组然后从目标分组里取出序列用多序列比对和系统发育推断还原它们的演化历史。在这个链路里OrthoFinder扮演的是“分拣员”角色它把来自不同物种的蛋白序列按直系同源关系聚成一个个Group简称OG保证你接下来建的树不会混入旁系同源基因而看不清真正的演化脉络。iTOL则是最后的“出版级修图师”负责把枯燥的newick文本变成带颜色、带分组、带支持值的可读图表。中间那几步——比对修剪和建树——才是决定基因树质量的关键很多人恰恰在这里图省事用了默认参数后就不管了。1.2 兰科案例的研究设定为了把整个流程讲具体我设计了一个可以直接套用的实际场景假设我们手头有5个兰科物种的蛋白组分别是小兰屿蝴蝶兰、深圳拟兰、铁皮石斛、天麻、白及再加两个外类群——芦笋和海枣。外类群的作用是给进化树定根并判断演化极性这点后面建树时还会细说。我们要研究的基因家族是MADS-box里的SEPALLATA类基因这类基因跟兰科高度特化的花结构特别是唇瓣和合蕊柱的发育有非常直接的关联。真实跑数据时你需要准备的输入就是每个物种一个蛋白序列fasta文件。文件命名建议全部用“属名_种名.fa”这种格式比如phalaenopsis_equestris.fa、dendrobium_catenatum.faID里只用英文字母、数字和下划线。OrthoFinder对输入文件名和序列ID里的特殊字符很敏感如果你图省事用了空格、括号、竖线甚至中文运行过程里往往不会立刻报错而是会产出一堆残缺结果排查起来非常痛苦。2. 数据准备与OrthoFinder运行手册2.1 输入文件的质量检查很多第一次跑OrthoFinder的朋友文件能读出来就直接开跑结果跑完之后发现某个物种的直系同源组数量异常少回头查才发现蛋白组里有大量低于50个氨基酸的碎片序列把同源推断结果搅得一团糟。我在做兰科数据时一般会先做两步检查。第一步是用SeqKit看每个文件的序列数量和长度分布seqkit stats phalaenopsis_equestris.fa dendrobium_catenatum.fa gastrodia_elata.fa第二步是过滤掉过短的序列。蛋白序列低于50个氨基酸基本不可能是完整功能蛋白保留它们只会干扰后续的Diamond比对和MCL聚类建议统一过滤掉seqkit seq -m 50 -g phalaenopsis_equestris.fa phalaenopsis_equestris.clean.fa-m 50表示只保留长度不小于50的序列-g会去掉序列里的终止密码子星号。过滤完最好再跑一次seqkit stats看一眼如果删掉的序列占比超过10%你得回头检查原始注释流程是不是出了问题。2.2 OrthoFinder的安装和核心参数安装OrthoFinder最省事的方式是用conda它会自动把依赖的Diamond、MCL、FastTree这些工具一起装好conda create -n orthofinder -c conda-forge -c bioconda orthofinder conda activate orthofinder用conda装的好处不只是省事它能保证OrthoFinder调用的Diamond版本和MCL版本是互相兼容的。我自己曾经图版本新手动装了最新的Diamond结果OrthoFinder的比对结果格式解析出现问题研究了半天才发现是版本匹配的锅。运行命令我常用这组orthofinder -f input_pep -t 20 -a 20 -M msa -o orthofinder_result参数说明如下-f指定存放所有物种蛋白fasta文件的目录-t是比对和聚类阶段用的线程数-a是物种树推断阶段用的线程数两者可以分开设置我一般直接都设成CPU核心数-M msa表示要做多序列比对这一步会顺带生成基因树如果你不加这个参数输出的结果里就没有Gene_Trees文件夹后续就没办法做基因树筛选了-o指定输出目录。一个完整的7个物种蛋白组跑下来每个物种大约3万到4万个蛋白20线程大概需要2到4个小时。运行过程中终端会打印当前阶段看到Clustering阶段结束后OrthoFinder会停下来在所有物种之间重新做一次配对推断这个过程不用干预等着就行。2.3 运行结果目录里到底有什么运行结束后OrthoFinder会生成一个类似Results_Nov12_2024_123456的目录里面最需要关注的是这几个Orthogroups/Orthogroups.tsv直系同源组总表行是OG编号列是物种单元格里是该物种在这个OG里的基因ID这是最核心的筛选入口。Orthogroup_Sequences/每个OG对应的原始蛋白序列文件文件名和OG编号一致后面建树直接从这里提取序列就行。Gene_Trees/每个OG的基因树newick文件文件名类似OG0016620_tree.txt。Species_Tree/SpeciesTree_rooted.txt物种树这是OrthoFinder依据所有单拷贝直系同源组推断出来的后面判断基因树拓扑是否合理时会用到。Comparative_Genomics_Statistics/各种统计表包括每个OG的基因数分布做大规模筛选时会用上。版本差异要特别注意OrthoFinder 2.x系列的输出和早期1.x系列差别很大比如早期版本里是Orthogroups.txt而现在主推的是带基因ID矩阵的Orthogroups.tsv另外现在2.5以上的版本还多了一个Phylogenetic_Hierarchical_Orthogroups/N0.tsv这是考虑了物种系统关系的层级正交组注释做深度的演化比较时会用到。如果你在网上下载的教程命令跟你本地输出对不上先确认OrthoFinder版本很多问题都是版本错位造成的。3. 从OrthoFinder结果中精准锁定目标基因树3.1 用Orthogroups.tsv做覆盖度和拷贝数初筛OrthoFinder默认会把所有物种共有的、结构清晰的单拷贝OG放在最前面但实际我们关心的基因家族往往不是单拷贝的。就拿SEPALLATA类MADS-box基因来说兰科基因组里通常有3到6个拷贝不同物种拷贝数还不一样这种有扩张的OG才是演化分析中最有意思的目标。筛选的第一步是快速扫描每个OG的物种覆盖度和总基因数。我一般用awk统计awk -F\t {n0; for(i2;iNF;i) if($i0) n; print $1, n, NF-1} Orthogroups.tsv og_coverage.txt sort -k2,2n -k3,3nr og_coverage.txt | grep -P \s7\s | head -n 50这段命令的意思是把第二列到最后一列里基因数大于0的物种数统计出来然后筛选那些覆盖了全部7个物种5个兰科2个外类群的OG。覆盖度全满是一个比较稳的筛选条件如果一个OG在某个物种里完全缺失你不清楚它到底是真实丢失还是注释漏了分析时很难解释除非你专门研究基因丢失。拷贝数方面不同物种的拷贝数如果差异过于悬殊比如某个OG里蝴蝶兰有8个拷贝而其他物种都只有1个这种强烈扩张提示可能存在串联重复或近期基因扩张事件很有价值但也意味着后续建树时要特别小心因为旁系同源基因混在一起基因树和物种树对不上是常有的事。3.2 看基因树初筛拓扑和功能注释双重确认初筛出候选OG后不要急着进下一步先把OrthoFinder生成的Gene_Trees里的基因树用任意能显示newick的工具打开瞄一眼判断这些基因在树上的聚类情况。你不需要分析得多精细至少确认目标OG里的基因不是完全随机堆在一起而是呈现出一定的谱系或物种聚类倾向。这一步能帮你提前发现一些明显有问题的OG比如序列嵌合体或者基因模型注释错误。如果条件允许再做一次功能注释做双重确认。OrthoFinder只管序列相似性和系统关系不负责告诉你这个OG是什么基因。我用的是eggNOG-mapperemapper.py -i OG0016620.fa --cpu 8 -o OG0016620_eggnog结果注释里如果明确出现MADS-box或SEPALLATA相关的条目这个OG就可以放心地进入建树流程了。如果你不想装eggNOG-mapper也可以用NCBI的CDD搜索网页版或者用本地InterProScan前者要手动一个个查后者跑得比较慢但都可行。3.3 提取目标OG的序列文件确认好目标OG之后直接从OrthoFinder的输出里提取序列是最省事的办法cp orthofinder_result/Results_Nov12_2024_123456/Orthogroup_Sequences/OG0016620.fa .如果你在Orthogroups.tsv里只能看到基因ID却找不到对应的序列文件那多半是你用了自己整理的矩阵来筛OG这时就得回到全蛋白组文件里用seqkit grep按ID批量提取cut -f1 OG_target_ids.txt | seqkit grep -f - phalaenopsis_equestris.clean.fa这里要注意序列IDS的一致性。Orthofinder在输出基因树时叶子节点标签保留的是输入fasta文件的序列ID所以你在iTOL里看到的名字多半是类似PH_eq0012345这样的ID而不是物种名加基因名的可读格式。为了避免后面iTOL里认不出物种我在进入建树前会把序列ID统一清洗成物种缩写|基因ID的格式用seqkit replace就能批量完成seqkit replace -p ^ -r Phal_equestris| phalaenopsis_equestris.clean.fa注意这里的竖线符号在后续iTOL注释文件里没问题但如果你要用一些命令行工具处理newick竖线偶尔会被当成特殊字符所以也有人用“”或者双下划线。我后来统一改用双下划线分隔最省心iTOL和大部分下游软件都不会误判。4. 多序列比对与IQ-TREE建树实操4.1 MAFFT比对别用默认参数直接跑拿到目标OG的序列之后第一件要做的事是多序列比对。这是整个流程里最容易被忽视、却最能决定后续建树质量的一步。很多人直接跑mafft input.fa output.aln就完事了但默认参数是兼顾速度和精度的折中方案对序列差异较大的跨物种同源基因来说并不够好。我跑基因家族分析时习惯用带迭代优化的L-INS-i策略它适合序列条数在200条以内的数据集准确度比默认的auto模式高不少mafft --localpair --maxiterate 1000 OG0016620.fa OG0016620.aln--localpair对应L-INS-i会做局部两两比对后逐步迭代优化--maxiterate 1000让迭代次数更充分。对于只有几十条序列的OG计算时间顶多几分钟。如果序列条数太多超过几百条L-INS-i会慢得让人失去耐心这时可以用--auto让MAFFT自己选策略。比对完成后别急着建树先花10秒钟检查比对质量。用alv或者Jalview打开比对结果重点看有没有大片的非保守区域和明显错位的gap。我见过太多人栽在这里序列比对里残留大片长度不一的non-conserved区域建出来的树上对应分支的支持值会异常低甚至整棵树的长枝吸引效应直接毁掉拓扑。4.2 trimAl修剪去掉“垃圾列”才能得到可靠树比对好的序列如果直接建树那些高度可变、gap率极高的区域会引入大量噪音信号相当于给进化分析加入了随机噪声。修剪这一步看着不起眼对树的支持值影响却很大。trimAl的-automated1参数会根据比对特征自动决定修剪阈值适合多数情况trimal -in OG0016620.aln -out OG0016620.trim.fa -automated1如果你做的是跨物种远缘比较序列差异比较大我建议再配合严格一点的gap和保守性过滤trimal -in OG0016620.aln -out OG0016620.trim.fa -gt 0.3 -st 0.001-gt 0.3要求该列至少有30%的序列不是gap-st 0.001要求该列至少保留0.1%的相似性这是一个平衡信息量和噪音的常用组合。修剪后的序列长度不要低于原始比对的50%如果修剪完只剩一小截序列说明原始比对质量出了问题应该回头改MAFFT参数而不是继续往下走。4.3 IQ-TREE建树模型选择和分支支持值一次说清建树工具我用IQ-TREE原因很简单它会自动用ModelFinder选择最适合的氨基酸替代模型而且内置了超快的自举法UFBoot2几千个位点跑得飞快。核心命令如下iqtree -s OG0016620.trim.fa -m MFP -bb 1000 -alrt 1000 -nt AUTO --prefix OG0016620逐项解释一下。-m MFP让IQ-TREE在众多替代模型里自动挑选考虑到兰科这些物种的序列相对保守选出来的模型很大概率是LG或JTT这类常用模型的变体-bb 1000是跑1000次UFBoot2快速自举这是评估分支支持度的主力-alrt 1000额外跑SH-aLRT检验两个支持率配合使用比单看一个可靠很多--prefix指定输出文件名前缀避免每次生成一堆unrooted这种默认名字。跑完以后输出文件里最需要关注的是.treefile和.iqtree。.iqtree是详细的报告里面有最终选的模型、位点信息、各分支支持值汇总我每次都会先扫一眼确认模型不是奇怪的混合模型.treefile是带支持值的newick树这就是你接下来要拖进iTOL的文件。判断支持率好坏有个通用标准SH-aLRT大于等于80且UFBoot大于等于95的分支算高支持两个都低的分支在描述结论时要特别谨慎不要拿着一个低支持分支去讲演化故事。这个标准不只适用这个案例做任何基因树我都建议照这个门槛来。5. iTOL可视化从上传树文件到发表级美化的完整操作5.1 树文件上传前的必要准备IQ-TREE生成的.treefile可以直接拖进iTOL但为了后续标注方便我会先确认一下叶子节点的ID格式。iTOL的注释文件是按叶子ID精确匹配的如果ID里有多余空格或者特殊字符注释文件会反复报错非常折磨人。我的建议是在建树之前就把ID规范好比如Phal_equestris_EQ12345这种前段是物种缩写后段是基因ID中间用双下划线连接。经过MAFFT、trimAl、IQ-TREE这一路ID不会变所以进iTOL之后就能直接用前缀区分物种。5.2 iTOL基础操作与注释文件格式打开iTOL网站建议先注册免费账号不然每次上传的树和注释文件都没法保存重新打开浏览器就要重新弄。上传时选择“Newick”格式把.treefile拖进去页面几秒钟内就会渲染出树。基础的样式调整在“Control panel”面板里枝的粗细调成2到3显示更清晰叶子标签字体大小调到20以上标签颜色可以按需设置。如果你想显示分支支持值在“Advanced”选项里勾选“Node labels”iTOL会默认读取newick里节点处的支持值。要注意IQ-TREE的.treefile里支持值写在节点位置的括号结构里iTOL读取时没问题但显示出来数字很小需要调大字体。分组上色是iTOL里最常用的操作。在Control panel的“Datasets”菜单里选“Import dataset”上传注释文件。我经常用两种数据集类型。第一种是DATASET_COLORSTRIP在树的外圈画彩色条带快速区分物种。格式如下DATASET_COLORSTRIP SEPARATOR TAB DATASET_LABEL Species COLOR #ff0000 COLOR_BRANCHES 1 STRIP_WIDTH 25 DATA Phal_equestris #ff0000 Phalaenopsis Den_catenatum #00aa00 Dendrobium Gas_elata #0000ff Gastrodia Blet_striata #ffaa00 Bletilla Apos_shenzhenica #aa00ff Apostasia Aspar_officinalis #888888 Outgroup Phoenix_dactylifera #444444 Outgroup这里的第二列是条带颜色第三列是图例里显示的文本。注意SEPARATOR TAB表示字段之间用Tab键分隔如果你从Excel里复制粘贴粘贴时极容易把Tab变成空格iTOL会直接提示解析错误。我吃过这个亏后来都是用文本编辑器写注释文件。第二种是DATASET_STYLE可以精确控制某个分支或单个叶子的颜色、粗细、甚至形状。比如我想把兰科的分支高亮成橙色DATASET_STYLE SEPARATOR TAB DATASET_LABEL Highlight COLOR #ff8000 DATA Phal_equestris_EQ12345 branch #ff8000 Den_catenatum_DC23456 branch #ff8000如果只想高亮一个clade而不是一条条列出所有叶子可以在iTOL页面里先用鼠标点选目标节点右键选择“Create dataset for selected branches”iTOL会自动生成包含该分支所有叶子的数据集你再微调颜色就行。5.3 树根位置与最终导出基因树默认是不定根的iTOL里展示时会自动取一个中间根。如果你的分析有明确的外类群比如这个案例里的芦笋和海枣应该把树根设在外类群的分支上。操作方法是在iTOL页面里点击外类群分支对应的节点右键选择“Root here”整棵树就会以该节点为根重排。这一操作对解读拓扑非常关键尤其是判断哪边是祖先状态、哪边是衍生状态。导出图片时我用得最多的是SVG和PDF矢量格式投稿期刊时放大多少都不糊。在“Export”菜单里选SVG分辨率选300dpi勾选“Scale bar”显示比例尺一件复杂的基因树图就完成了。如果你要在后续PPT里继续编辑SVG可以直接拖进Illustrator或者Inkscape改细节。iTOL还有一个容易忽略的实用功能——折叠大分支。如果某个clade特别大挡住了焦点区域右键选中该clade的根节点选“Collapse”整枝会收成一个三角树图瞬间清爽很多。这在展示超大基因家族树时几乎是必需操作。6. 常见问题与排查技巧实录6.1 OrthoFinder运行阶段的问题OrthoFinder跑挂的最常见原因是内存不足。7个物种、每物种3万条蛋白实际内存大概需要20到30GB如果你是在8GB内存的笔记本上跑我建议缩小输入规模比如先用较大的isoform代表去冗余或者减少物种数。另外确认-t线程数不要超过物理核心数线程设得过多反而会让Diamond的比对阶段频繁分配内存OOM报错概率更高。还有一个隐蔽问题输入目录里的隐藏文件。OrthoFinder会扫描-f目录下的所有文件如果目录里有.DS_Store、*.txt说明文件等杂项它会尝试当fasta解析然后报错。解决办法是输入目录里只放每个物种一个.fa文件不要放其它任何东西。我习惯把原始数据另外备份输入目录保持干净。6.2 筛选OG和建树阶段的问题很多人在筛选OG时只看某一物种的拷贝数忽略了全物种覆盖度结果选中一个在模式物种里有3个拷贝、但外类群缺失的OG建出来的树缺少外群定根只能靠midpoint解释结果时底气不足。所以筛选时一定以“全物种覆盖优先拷贝数差异其次”为原则。建树后出现基因树与物种树明显不一致是最让新手恐慌的情况。比如天麻的序列没有跟天麻物种所在分支聚在一起而是跟蝴蝶兰聚在一起。先别急着怀疑流程错了这一般说明该OG里包含旁系同源基因也就是不同拷贝的分化发生在物种分化之前基因树里呈现的是基因重复事件后的演化关系。这时你有两个选择将不同拷贝拆开重新做树或者直接以基因树为准来解读基因扩张事件。只要支持值可靠基因树与物种树的这种差异本身就是重要的演化证据。6.3 iTOL注释文件报错速查iTOL注释文件报错的九成原因是分隔符问题。记住SEPARATOR TAB之后DATA区每一行的列之间必须是真正的Tab字符。有些编辑器看上去是Tab实际是连续空格Windows记事本编辑过的文件还可能带\r回车符。建议统一用VS Code或者Sublime Text把文件编码保存为UTF-8 without BOM。如果注释文件里写的是叶子ID不匹配iTOL会在该行报错“Cannot find node”。这多半是因为你改过序列ID但忘了同步更新树文件或者建树过程里某个工具悄悄在ID里加了后缀。排查办法是把.treefile下载下来搜一下ID和注释文件逐字对照肉眼通常几秒钟就能发现问题。另一个容易踩的点是iTOL数据集类型名字不能拼错DATASET_COLORSTRIP少写一个字母都会直接报错复制模板改比手打靠谱。6.4 支持率显示与图形导出的几个坑IQ-TREE生成的.treefile里支持率都在但iTOL默认不显示节点标签很多人以为支持率丢了。在Control panel里找到“Node labels”选择“display”支持率就出来了。如果你的树特别大节点支持率密密麻麻可以在“Advanced”里调整只显示大于某个阈值的支持率比如只显示UFBoot大于等于95的节点图更干净也更有说服力。导出图片时如果发现叶子标签被截断多半是画布尺寸没调好。iTOL的导出设置里可以手动输入画布宽度或者勾选“Fit to page”。另外SVG目录下如果文件巨大可以先在iTOL里折叠那些不重要的clade再导出文件体积会小很多处理起来也更流畅。这套流程我前前后后跑过不下十遍踩坑最多、最不值得的永远是两类一类是不看输入文件直接开跑另一类是序列ID格式不一致导致后期处处碰壁。OrthoFinder和iTOL都是非常成熟稳定的工具只要按部就班把输入规范好把每一步的参数意义理解到位整个流程跑下来非常顺。最后再分享一个小技巧完成一次分析后把所有的命令、版本号、筛选条件原样保存到一个Markdown文件里哪怕只是自己看也能在半年后写论文方法部分时替你省下大量回忆时间。环境是会变的数据是会重跑的但记录下来的命令和方法不会背叛你。