做基因家族分析的时候经常会有人问我系统发育树画出来了基因结构也展示了怎么证明这些基因在不同物种之间是真的同源而不是我强行把它们拉到一起我的做法是补一张物种间共线性分析图。共线性分析synteny/collinearity analysis能说明两边的同源基因不仅序列相似而且它们周边的基因排列顺序也保持保守这个证据比单独一条BLAST结果要硬得多。最早做这件事一般都要在Linux下装MCScanX再用Perl脚本清洗数据对没接触过命令行的朋友不太友好。后来我长期用TBtools它把MCScanX这一整套流程封装成了图形界面双击就能跑还能直接出图。这篇文章就按我实际跑的完整流程从输入数据准备到最终出版级图片一步步拆开说清楚适合正在做基因家族、比较基因组和进化分析的朋友参考。1. 物种间共线性分析到底在回答什么科学问题1.1 共线性不是简单的“序列相似”很多初学者会把共线性理解成“两个物种的序列很像”其实不对。序列相似是BLAST就能回答的问题共线性回答的是另一个问题这些相似的基因在两条染色体上的相对顺序是不是也一致。我常用一个生活化的类比BLAST相当于告诉你两栋楼里都有“客厅”和“卧室”共线性则进一步说明客厅旁边是卧室卧室旁边是厨房两栋楼的房间布局都差不多。如果两个物种在物种分化之后没有发生大规模染色体重排那么它们从祖先基因组里继承下来的那些基因通常会以“模块”的形式保留在染色体上顺序大致相同。这个模块就是共线性区块syntenic block / collinear block。在比较基因组学里synteny更强调区块层面的同源关系collinearity则更强调区块内基因顺序的保守。不过日常交流中大家经常混用TBtools里叫“共线性分析”实际上跑的是MCScanX的collinearity检测。1.2 物种间共线性与物种内共线性的用途差异我经常被问到物种间共线性和物种内共线性有什么区别这两个功能在TBtools里都能做但回答的科学问题完全不同。物种内共线性intra-species collinearity用于检测单个基因组内部的重复区块典型场景是全基因组复制事件WGD分析。比如拟南芥经历过多轮WGD你拿它自己的染色体和自身比对能看到很多同源区块分布在不同的染色体上这些区块就是古代加倍事件留下的痕迹。物种间共线性inter-species collinearity用于比较两个或更多物种之间的保守区块。做基因家族分析时如果某个家族成员在物种A和物种B里位于彼此对应的共线性区块中就可以很有底气地说这两个基因是直系同源至少也是在物种分化之前就从共同祖先继承下来的而不是各自独立扩张出来的。这个证据对后续的基因命名、功能注释和进化推断都很有帮助。1.3 一张共线图能读出多少进化信息完成一次物种间共线性分析后你能从结果里看到的东西比想象中多哪些染色体区段在物种分化后保持完整哪些被倒位、易位打断哪些区块里的基因发生了丢失甚至可以根据同义替代率Ks估算大致的分化时间。这些信息在系统发育树之外提供了一个独立的证据维度。举个例子水稻和高粱同属禾本科亲缘关系较近它们的共线性区块应该非常完整一条水稻染色体可能对应一到几条高粱染色体。而水稻和拟南芥这类单双子叶之间的比较共线性区块会比较破碎需要更宽松的阈值才能看到残留的保守片段。所以拿到结果后不要急着下结论先看看两个物种的亲缘关系再判断结果是否符合预期。2. 动手前先分清TBtools里的分析工具和绘图工具2.1 One Step MCScanX 才是真正做共线性分析的工具TBtools里和共线性有关的功能很多我见过不少新手把Dual Systeny Plot当成“共线性分析”直接点开然后发现不知道怎么填文件。这里先明确一下分工真正负责计算共线性区块的是One Step MCScanX其他多数工具都只是对MCScanX结果进行可视化。One Step MCScanX顾名思义一步完成“蛋白序列BLAST比对 MCScanX共线性搜索”。MCScanX本身是一个命令行程序原理是基于BLAST得到的同源基因对再用动态规划算法把那些在染色体上顺序一致、方向一致的同源基因对连成区块。对于不熟悉Linux的人来说光是把输入文件整理成MCScanX要求的格式就够折腾的TBtools在这里的价值就是把格式处理、程序调用和结果整理全部包办了。我自己的习惯是只有在需要批量处理几十个物种时才会回到命令行去跑MCScanX原版。普通的两两比较直接在TBtools里点一下就完成了保证结果一致还省去环境配置的麻烦。2.2 可视化工具Dual Systeny Plot和Advanced Circos的分工MCScanX算完之后得到的是一堆包含“染色体A上的基因列表 vs 染色体B上的基因列表”的文本结果直接看很难形成直觉所以要画图。TBtools里常用的有两个画图工具。工具适用场景输入内容输出效果Dual Systeny Plot两个物种之间的共线性展示两个物种的GFF文件 MCScanX的collinearity结果上下两条染色体/图谱之间用曲线连接共线基因对Advanced Circos多物种、多区块的环形展示染色体长度文件 共线性block文件 可选基因标注环形圈图外圈染色体内部连线表示共线区块从项目沟通角度来说给合作者或审稿人看时圈图Circos往往最直观因为可以把多个物种、多个共线区块放在一张图里信息密度很高。但如果只需要简单看看两个物种之间的关系Dual Systeny Plot更轻量出图也更快不用花时间调整一堆圈图参数。2.3 我建议的工作流我在实际项目中的工作流一般是四步第一步先做一个小规模测试。比如只取目标物种的某一条染色体和另一个物种做One Step MCScanX或者用完整数据但设置较严格的E-value先跑一次。这么做的目的是快速验证GFF和蛋白序列的ID格式有没有问题避免等了两小时后才发现文件前缀不匹配。第二步跑全基因组的One Step MCScanX得到正式的collinearity结果。第三步根据需求画图。两个物种用Dual Systeny Plot超过两个物种或者要给文章做摘要图就用Advanced Circos。第四步从collinearity文件里提取具体基因对结合基因家族成员列表做筛选把特定家族的共线关系单独提出来再叠加标注到图上。别看这四步简单每一步都有很多细节尤其是第一步的数据准备最不起眼也最容易翻车。3. 输入文件准备GFF和蛋白序列的规范化处理3.1 从公共数据库下载哪些文件做物种间共线性分析需要每个物种提供两个核心文件一个是基因组注释文件GFF/GFF3另一个是蛋白序列文件FASTA。GFF文件里记录的是每个基因在染色体上的位置、方向、转录本结构等信息。蛋白序列文件则提供了所有基因对应的氨基酸序列。下载来源我一般首选Ensembl Plants、Phytozome或NCBI RefSeq这几个数据库的物种覆盖面广注释质量也比较可靠。如果你做基因家族分析前面可能已经下载过目标物种的蛋白序列那这时候只需要再补一个GFF文件。注意下载时看清是GFF3还是GTF推荐GFF3。GTF主要给转录本定量用缺少MCScanX需要的部分层级信息即使能跑也容易出问题。下载后先解压然后用文本编辑器或less命令瞄一眼文件内容确认没有空文件、没有乱码再继续下一步。3.2 文件命名和目录规范这一步我吃了不少亏现在会特别强调TBtools对中文路径和文件名中的特殊字符非常敏感。建议创建一个纯英文目录例如D:/synteny/把两个物种的文件放进去并统一使用简洁前缀命名。我常用的命名规则是“物种拉丁名前几个字母缩写 文件类型”例如水稻Osativa.gff3、Osativa.pep.fa高粱Sbicolor.gff3、Sbicolor.pep.fa其中Os是Oryza sativa的缩写Sb是Sorghum bicolor的缩写。文件名里只保留字母、数字、点和下划线不要有空格也不要以数字开头。原因很简单MCScanX在运行过程中会用文件名前缀作为样本标识如果前缀很乱后面解析结果时你会非常痛苦。提示文件名前缀在同一轮分析里不要重复。比如同时做水稻和水稻的物种内共线性分析两个文件不能都叫rice.gff3否则输出文件会互相覆盖。用不同的拉丁名缩写最稳妥。3.3 GFF与FASTA的ID对应关系检查这是整个流程里最关键的检查点没有之一。One Step MCScanX需要把蛋白序列和基因位置对应起来对应关系就是FASTA序列头的ID和GFF文件中基因/转录本的ID。不同数据库的ID风格差异很大。Ensembl Plants通常用gene:...或转录本ID作为FASTA序列头NCBI则经常出现lcl|前缀。如果GFF中的基因ID是Os01g0100100而FASTA头的ID是XP_015612345.1那MCScanX就会认为它们毫无关系最后结果基本是空的。我检查ID是否匹配的办法很简单用文本编辑器分别打开GFF和FASTA看两种ID的写法是否一致。或者用TBtools里的GFF3 to Table工具把GFF转换成表格和FASTA头部做个对比。3.4 我常用的预处理命令如果发现ID不匹配最简单的处理是通过命令行批量改名。下面这两个命令我经常用# 查看FASTA前10条序列的ID grep ^ Osativa.pep.fa | head -n 10 # 查看GFF中基因特征的前5行 grep -P \tgene\t Osativa.gff3 | head -n 5如果FASTA头里除了ID还有其他描述信息需要只保留第一个字段# 将序列头统一改成第一个下划线/空格前的ID例如gene_001 description - gene_001 awk {print $1} Osativa.pep.fa Osativa.clean.pep.fa注意awk {print $1}这种写法对标准FASTA头有效对NCBI那种lcl|xxx开头的头也能把描述去掉。如果要去掉lcl|前缀可以配合sed一起用sed s/^lcl|// Osativa.ncbi.pep.fa Osativa.clean.pep.fa改完之后一定要再看一眼确认ID之间没有多余的空格或特殊符号。这个环节多花十分钟后面能少折腾两个小时。4. 核心一步One Step MCScanX参数逐项讲解与运行4.1 工具入口和界面TBtools启动后在顶部菜单栏找到Comparative Genomics里面通常有One Step MCScanX不同版本的位置可能略有差异有的放在Genomics菜单下但名称不变。点击后会弹出一个窗口需要填写的核心内容有三个GFF文件、蛋白序列文件和输出目录。窗口里的文件选择顺序是先选参考物种的GFF和蛋白序列再选第二个物种的GFF和蛋白序列。严格来说MCScanX的运行不区分“参考”和“查询”它会把两个物种的蛋白放在一起做双向BLAST再搜索共线性区块所以这里的顺序只影响输出文件的组织结构不影响最终是否发现共线基因对。输出目录建议单独建一个文件夹比如./synteny_out不要直接输出到桌面上也不要用中文路径。TBtools会把BLAST中间文件、MCScanX结果文件都写到这个目录里目录越干净越容易整理。4.2 核心参数E-value和Max GapsOne Step MCScanX界面里有两个参数我每次都会调整E-value和Max Gaps。E-value是BLAST筛选同源基因对的阈值默认通常是1e-5。E-value越小筛选越严格得到的同源基因对越少但越可靠E-value越宽松越能找回远端同源基因但也会引入更多假阳性。我在实际操作中的经验是近缘物种比如禾本科内部用默认的1e-5就很好远缘物种比如动物和植物之间可以放宽到1e-10或1e-5之间重点看是否能找回已知的保守基因如果结果里共线性区块太多太碎可以收紧到1e-7试试。Max Gaps参数很容易被忽略它指的是在共线性区块中两个相邻同源基因之间允许插入多少个无关基因。MCScanX的经典默认阈值是25。可以这样理解如果两个物种在某个区段的共线性基因排列得非常密集中间没有太多其他基因说明这个区块非常保守如果中间塞了几十个不相关基因那这段“共线性”可能只是偶然的相似可信度不高。MCScanX还有一个概念叫anchors锚定基因也就是处在共线性区块里的同源基因对。使用默认参数时一个共线性区块至少要包含5对anchor基因才会被记录。这个最小anchor数量在部分TBtools版本中也可以调整如果两个物种分化太远可以适当降到3但我不建议降到2以下否则很容易从随机BLAST结果里拼出大量假区块。4.3 输出目录与运行监控点击运行后TBtools底部会弹出日志窗口显示当前进度。One Step MCScanX内部做的事情大致如下解析两个物种的GFF文件整理出每个基因的染色体位置调用BLASTP对两个物种的蛋白序列做两两比对将BLAST结果转换成MCScanX要求的格式运行MCScanX搜索共线性区块整理输出文件。日志里如果出现blastp not found之类的提示说明系统里没有安装NCBI BLAST或者TBtools没有正确找到BLAST程序。解决办法是提前安装BLAST并把blastp所在目录加入环境变量PATH。在Windows下如果你用Anaconda也可以打开Anaconda Prompt执行conda install -c bioconda blast然后重启TBtools让TBtools继承这个环境变量。运行时间取决于物种基因组大小和蛋白数量。两个蛋白组各3万个基因的物种完整跑完通常需要几十分钟到数小时不等主要瓶颈在BLASTP这一步。所以我在第二节提到的小规模测试非常重要先用小数据把流程跑通再上全基因组能帮你省下不少等待时间。4.4 输出文件清单运行结束后输出目录里会出现一系列文件。其中最重要的两个.blast或类似命名的BLAST结果文件记录所有同源基因对.collinearity文件记录最终检测到的共线性区块。另外还会有一些MCScanX产生的辅助文件和TBtools的日志文件。我一般不会直接删掉因为后续统计、画图常常会用到。拿到.collinearity文件后可以先用文本编辑器打开看看确认里面不是空的再做可视化。5. 出图环节从Dual Systeny Plot到Advanced Circos5.1 Dual Systeny Plot两物种点图怎么调参数当我只需要快速查看两个物种之间的共线性情况时首选Dual Systeny Plot。它画出来的效果是上方一条参考物种的染色体图谱下方一条查询物种的染色体图谱中间用曲线把共线基因对连起来。打开TBtools后在Comparative Genomics菜单下找到Dual Systeny Plot。界面里需要提供参考物种的GFF文件查询物种的GFF文件MCScanX得到的.collinearity文件。填完之后点击StartTBtools会根据GFF里的染色体长度和collinearity文件里的基因对位置自动生成一条竖直方向的连线图。默认情况下每条染色体用不同颜色区分共线基因对之间用灰色曲线连接。如果图里的连线过于密集可以调整一个显示比例参数只显示那些anchor数量较多的共线区块把碎片化的单基因连线过滤掉。这个阈值设置在界面上通常叫类似“min score”或通过滑块控制具体名称随版本不同思路是把可信度低的支行去掉让主干结构更清晰。从图里能直观看到三种信息一是共线区块所在的染色体区域二是倒位事件区块方向反转连线呈X形交叉三是转座事件同一段序列连接到了另一条染色体的不同位置。5.2 Advanced Circos多物种环形图的配置逻辑如果要展示三个以上物种或者想在图上同时标注基因家族成员我会用Advanced Circos。这个工具功能非常强但配置比Dual Systeny Plot复杂一些因为它本质上是一个自定义程度很高的圈图绘制器。基本的配置步骤是这样在Graphics菜单下找到Advanced Circos然后在界面里分几个区域添加数据文件第一个区域是染色体长度文件。格式是两列第一列染色体ID第二列染色体长度用Tab分隔不要加表头。这个文件可以从GFF里提取也可以自己在Excel里整理后另存为txt。第二个区域是共线性区块文件。可以来源于MCScanX的.collinearity也可以使用BLAST结果转换而来的区块文件。如果选.collinearityTBtools通常能直接识别并解析如果同时展示多个物种之间的共线性需要把多组区块文件都导入进来。第三个区域是可选的标注文件用来在圈图上高亮特定基因。做基因家族分析时这一步会非常加分把你关注的家基因成员ID整理成一个列表导入后在圈图对应位置就会标出名字或符号审稿人一眼就能看到该家族哪些基因位于共线区块中。5.3 出图与导出不管用Dual Systeny Plot还是Advanced Circos最后导出图片时我都会设置一个较高的分辨率。TBtools的绘图窗口通常支持调整宽度、高度和字体大小可以在界面里先把颜色、线条粗细调到满意再一次性导出PNG或PDF。建议有条件的话导PDF后续用AI或Inkscape做细节修改不会糊。另外不要忽略配色。一张好的共线性图信息量再大配色如果太乱也会显得很不专业。TBtools的默认配色已经经过设计日常够用如果要投稿我会把染色体颜色改成同一色系的深浅渐变共线连线的颜色改成半透明的灰色这样整体看起来更清爽。6. collinearity结果文件解读与后续常用分析6.1 .collinearity 文件到底长什么样MCScanX的.collinearity文件格式有固定结构打开后大概是这个样子## Alignment 0: score1254.0 e_value0.0 N12 L0 C0 0 0: gene_a_1 gene_a_2 gene_a_3 gene_a_4 ... 0 1: gene_b_1 gene_b_2 gene_b_3 gene_b_4 ... ## Alignment 1: score986.0 e_value0.0 N10 L0 C0 ...头一行里N后面的数字表示这个共线性区块里含有的基因对数量。下面两行一行来自参考物种一行来自查询物种每一对对位位置的基因就是一对共线基因。方向由和-表示如果两行方向一致说明区块没有发生倒位如果方向相反说明其中一个物种里发生了倒位。我以前刚接触这个文件时习惯直接用Excel打开但大文件很容易卡死后面改用文本编辑器查看配合命令行统计效率高很多。6.2 如何统计共线基因对数量统计共线基因对之前先要明确统计口径是统计所有区块里的总基因对数还是只统计某个基因家族成员涉及的共线关系。后者在基因家族分析中更常用。假设你的目标家族成员ID清单在family_ids.txt里可以用一个简单的循环来统计# 从collinearity文件中提取所有基因ID grep -v ^# Osativa_Sbicolor.collinearity | \ awk {$1; print} | tr \n | \ grep -v ^$ | sort -u all_genes_in_collinearity.txt # 和家族成员清单取交集 comm -12 (sort family_ids.txt) (sort all_genes_in_collinearity.txt) family_genes_in_collinearity.txt拿到交集列表后再回到原始collinearity文件里检索每个家族基因位于哪个区块、和另一个物种的哪个基因配对。这一步用Excel的筛选功能也能做但基因对一多就容易乱我一般会写个简单脚本处理。6.3 结合Ks值判断复制事件时间共线性区块只告诉你基因顺序保守没告诉你这些复制事件发生在多久以前。想要估计时间常用的指标是Ks同义位点替代率。Ks越小说明两个基因分化时间越近Ks越大说明分化时间越古老。TBtools里有个Simple KaKs Calculator可以批量计算一对基因的Ka和Ks。流程是先提取所有共线基因对的CDS序列按“reference_gene和query_gene”的配对关系整理成输入文件然后用Simple KaKs Calculator批量计算。算完之后把Ks值分布做成柱状图就能看到是否存在一个比较集中的近期扩张峰还是整体都是古老事件。对于做基因家族分析来说这个分布直接关系到你对“家族扩张机制”的判断——如果某个家族成员大量出现在Ks值较小的共线区块说明这个基因家族在近期发生过扩张可能和物种适应性进化有关。6.4 从物种间共线性走进基因家族分析提到基因家族分析是因为物种间共线性分析实在是一个高频前置步骤。常见的情形是你已经鉴定出一个目标基因家族在拟南芥里有10个成员在水稻里有15个成员接下来需要回答水稻多出来的5个成员是物种分化后自己扩张出来的还是拟南芥丢失了这时物种间共线性分析就能帮上忙。把目标家族的成员ID和.collinearity文件中的共线基因对做交集如果水稻某个成员在拟南芥中有一一对应的共线基因说明它在物种分化前就存在如果水稻有几个成员位于同一个共线区块中而拟南芥对应位置只有1个基因那就说明水稻这个分支上发生过串联复制或片段复制家族的扩张机制就有了比较坚实的证据。这也是为什么很多基因家族方向的文章里即使不专门讲比较基因组也会放一张共线性圈图的原因——它是逻辑链上非常有力的证据。7. 典型踩坑实录与排查思路7.1 报错“File not found”或“Unsupported GFF format”这个报错我遇到太多次了。绝大多数时候不是文件真的不存在而是路径里有中文、空格或特殊字符TBtools在解析时找不到正确位置。解决办法很简单把整个分析目录移动到C:/或D:/下一个纯英文路径文件名也改成英文重新运行。另一个常见原因是GFF文件本身已经不是标准格式。有些数据库下载的GFF注释文件行内用了空格而不是Tab做分隔符或者坐标列的顺序不对。TBtools会直接拒绝。检查方式是用文本编辑器打开GFF确认每一列之间是Tab最后一步再用GFF3 to Table工具转成表格看看是否正常。7.2 运行结束了但collinearity文件为空运行结束没有报错结果文件却是空的优先级最高的排查项就是GFF和FASTA的ID不一致。MCScanX是靠ID把蛋白序列映射回染色体位置的只要有一批基因的ID对不上这些基因就会被丢弃剩下能参与分析的基因越来越少自然检测不到区块。我的排查步骤很固定打开.collinearity文件看是否完全为空还是只有几个区块打开GFF文件随机挑几个基因ID打开FASTA文件搜索同样的ID如果不一致用awk或sed统一修改FASTA头部。另外还有一种情况两个物种亲缘关系特别远同源基因的BLAST E-value全部大于默认阈值导致找不到可用同源基因。这时把E-value放宽到1e-3可以救回来一部分但随之而来的假阳性也要警惕。7.3 结果太少不要急着怀疑软件很多人跑完发现共线性区块很少第一反应是TBtools出问题了。实际上区块少往往说明两个物种之间确实经历过较频繁的重排或者它们的亲缘关系确实较远。这时候先做两件事第一从BLAST结果文件里统计一下两个物种到底有多少对同源基因。如果同源基因对本身就很少那共线性区块少是正常结果调整参数意义不大。如果同源基因对很多但共线性区块很少说明这些同源基因分散在染色体各处没有形成保守顺序这时候可以把Min Anchors从5降到3看能不能恢复一些短区块。第二检查一下两个物种的染色体命名方式是否一致。有时候数据库里一个物种用chr1另一个用1MCScanX虽然能处理但画图时会出现染色体名对不上的问题看起来像结果里少了一片。7.4 绘图时染色体名对不上Dual Systeny Plot和Advanced Circos对染色体名的匹配非常严格。比如GFF里写的是Chr01collinearity文件里却写成chr1连线就会丢失。遇到这类问题建议从头统一染色体命名在GFF中用查找替换把Chr01改成chr1同时重新跑一遍One Step MCScanX因为collinearity文件里的染色体名是从GFF解析来的只改可视化文件不动源头问题仍会反复。7.5 大基因组跑不动先做“代表转录本”不少植物基因组动辄5万个基因而一个基因往往有多个转录本如果蛋白FASTA文件里把所有转录本的蛋白都放进去BLASTP的比对次数会急剧膨胀运行时间会变成“天”级别。我的做法是先提取每个基因的“代表转录本”也就是最长转录本或最长蛋白序列。这个操作在TBtools里有现成工具也可以自己写脚本。提取最长转录本的原则其实很简单按基因ID对转录本分组保留CDS或蛋白序列最长的一条然后重新生成FASTA文件。GFF这边也需要同步过滤只保留代表转录本对应的注释行。如果你觉得处理转录本麻烦还有一个替代方案在Ensembl Plants下载时直接选择“pep.all.fa”里面包含了所有转录本另外有些网站提供“representative transcript”或“primary transcript”版本的蛋白文件优先选这种。提示跑全基因组之前务必先拿一条染色体做测试。染色体级别的数据几分钟就能出结果既能验证ID、路径、BLAST环境都没问题又能让你提前看到collinearity文件的大致格式正式跑的时候心里有底。最后再说一点实际操作中的体会用TBtools做物种间共线性分析真正的门槛从来不是点击那几下鼠标而是进入软件之前的数据准备和离开软件之后的结果判断。我在实际项目中养成的习惯是每次下载数据后先把GFF和FASTA的ID对应关系检查一遍再顺手给文件改好规范命名最后用一条染色体做小规模测试。这套流程看起来多花了大半个小时但换来的是一整天不折腾。另一个小技巧是保存好中间文件。很多人跑完One Step MCScanX就把BLAST结果和日志删了只留一个collinearity文件后面想检查某个基因对是怎么比出来的只能重新跑一遍很浪费时间。我的习惯是把每个步骤的输入输出都放在同一个项目文件夹里按01_rawdata、02_input、03_result、04_figure这样分目录排好之后写文章、补实验、回应用审稿人意见时都能在几分钟内找到对应文件。最后想说共线性分析的结果不要孤零零地看。它最好配合系统发育树、基因结构、表达谱和Ks值一起使用才能把一个基因家族的进化故事讲完整。希望这篇东西能帮你少踩一些我踩过的坑顺利产出你自己满意的图。