每年都会有做群落生态学的朋友被同一个问题卡住样地调查的物种名录早就整理成了Excel接下来要算系统发育多样性PD、净亲缘关系指数NRI、系统发育信号这些指标可是系统发育树从哪来实验室没有测序条件也不可能一个种一个种去NCBI下载序列、比对、建树。我师弟去年年底就抱着300多个种的样地名录来找我愁眉苦脸地问怎么才能快速拿到一棵能用的树。我给他的方案很直接用V.PhyloMaker从物种列表直接生成一棵覆盖全部目标物种的维管束植物系统发育树全程只需要一份格式规范的物种名录在R里跑几分钟就出结果。V.PhyloMaker是R环境里的一个系统发育工具包方法论文发表在Methods in Ecology and Evolution上Jin Qian, 2022。它做的事情一句话就能说清楚把用户提供的物种学名列表通过名称匹配的方式挂载到一个已经构建好的、包含数万种维管植物的大尺度时间校准系统发育树通常称为mega-tree超树上然后抽出一棵只包含目标物种、且带有分支长度的系统发育树。如果你正在做群落构建、生物地理格局、性状进化这类研究手里有物种清单但缺少分子数据这篇文章应该能帮你省下大量的时间。1. 先搞清楚V.PhyloMaker到底在做什么1.1 不是所有系统发育树都需要自己建常规的建树流程是什么样的相信做过的人都知道先从每个物种上获取DNA条形码或者叶绿体基因组片段做多序列比对然后确定最适合的核苷酸替代模型再用最大似然法或贝叶斯推断重建拓扑结构还要做化石校准或者地质事件校准来获得分化时间。这一套流程走完顺利的话一两周不顺利的话光是序列获取和比对质量就能把人磨到没脾气。而且很多做生态学研究的人并不需要自己从头去推断演化关系——过去二十多年已经有大量高水平的系统发育研究把维管束植物的主干拓扑和主要支系的分化时间梳理得相当清楚。V.PhyloMaker的思路就是把这类公开成果借过来用。你不需要提供任何分子序列只需要提供一份物种名录。包里内置了一个巨大的、已经带有分支长度的超树覆盖了几万种维管植物。你的任务就是让物种名录里的每一个名字在这棵超树上找到自己该在的位置。这个过程在生态学里非常实用尤其是做区域尺度的群落系统发育结构分析时我们需要的往往不是某个基因片段推断出来的精细拓扑而是一个建立在学界共识之上的、可复现的演化框架。V.PhyloMaker干的就是这件事。1.2 超树加嫁接一个朴素但很有效的思路超树可以理解成一个巨大无比的家谱里面记录了数万种植物之间的亲缘关系以及大致分化时间。V.PhyloMaker拿到你的物种名录后会逐一做名称匹配如果某个物种本身就在超树里直接保留如果不在就根据它所属的属和科把它嫁接到超树对应的分支上。这个嫁接的位置通常是该属或科的基部节点附近相当于在属的旁边增设一个新分支。分支长度则是从超树继承的保留了演化时间的尺度所以最终拿到的树带有一致的时间信息。这个机制也决定了它的适用边界它不追求基于你样本自身序列证据的精确推断而是提供一个基于当前最大规模学术共识的近似演化框架。对绝大多数宏观生态学和生物地理学问题来说这个尺度的树已经足够用了。我自己用过好几次只要物种名录干净、分类系统统一生成的结果和后续用分子序列单独构建的核心树在属级以上拓扑上几乎一致。当然如果研究对象集中在某个物种极其丰富、且存在大量近期分化的类群里那就需要谨慎一点超树的分辨率可能不够这时候还是应该走传统建树流程。1.3 V.PhyloMaker、PhyloMaker、V.PhyloMaker2怎么选很多人第一次听到这个名字会有点晕因为相关工具不止一个。早期有个PhyloMaker2015年用的是当时版本的维管植物巨树但默认树已经比较旧了。后来出现了基于更新、更大规模超树的V.PhyloMaker也就是本文主要使用的版本。再后来又有V.PhyloMaker2额外支持用户提供自己的构建树比如你自己建好的某科系统发育树把自定义树作为嫁接的骨架。对绝大多数第一次使用的人我建议直接装V.PhyloMaker因为它的内置超树已经很大覆盖范围足够广用法也最简单。等以后你有了自己特定的关注类群、想把手头已有的某棵子树整合进来时再考虑V.PhyloMaker2也不迟。2. 开工前的准备R环境和物种名录格式2.1 安装包和基础依赖V.PhyloMaker的运行环境很常规普通R环境就能跑。安装命令很直接install.packages(V.PhyloMaker) library(V.PhyloMaker)如果CRAN安装失败最常见的原因是R版本偏旧可以先升级R到4.x再试。实在不行就通过GitHub安装作者仓库里的开发版。包本身不大但运行时会加载内置超树数据内存占用会上去我建议电脑内存至少8GB加载期间不要同时开一堆大型软件否则可能会卡死。第一次运行的时候我习惯先执行一下library(V.PhyloMaker)看看有没有报错再检查一下内置数据是否存在data(package V.PhyloMaker)正常情况下应该能看到类似GBOTB.extended的超树数据名称。这一步花不了几秒但能提前发现依赖缺失的问题比后面中途报错再回来排查要省心得多。2.2 spList三列格式species、genus、familyV.PhyloMaker的输入是一个数据框最常见的形式是三列species、genus、family。第一列是完整的学名也就是属名种加词比如Fagus longipetiolata第二列是属名第三列是科名。为什么科名必须给因为名称匹配是分优先级进行的物种名如果能精确匹配上超树直接保留匹配不上就用第三优先级里的科名来兜底把该物种挂到科内合适的位置。这份表是后续所有操作的基础格式越规范成功率越高。这里有个经常被忽略的细节species列必须是两段式学名格式不要带sp.、cf.这些限定词也不要把命名人比如Linn.写进去。V.PhyloMaker对种下等级的支持非常有限如果你研究的是亚种、变种建议统一合并到种一级来分析否则很可能会出现大量匹配失败。2.3 学名清理哪些坑必须在录入之前填平学名清理是整个流程里最花时间的一步也是决定匹配率的胜负手。我自己踩过的坑主要有三类。第一类是异名问题。分类学家隔几年就会修订一次系统发育框架以前叫Cyclobalanopsis glauca的现在很多系统里已经并入Quercus以前叫Michelia的一些种现在被并入Magnolia。如果名录里仍然使用旧学名超树可能根本找不到对应的物种节点。处理办法是找一个权威植物名录去做有效名校对比如POWOPlants of the World Online、中国植物志的在线系统或者直接用R包taxize批量校对。分类系统统一之后匹配率通常会明显提升。第二类是科名的系统归属问题。过去教科书上把枫香树放在金缕梅科Hamamelidaceae但APG IV系统已经把它放入枫香科Altingiaceae过去把槭树属单独放在槭树科Aceraceae现在则归入无患子科Sapindaceae。V.PhyloMaker内置超树使用的是较新的APG框架所以我在准备名录时会统一按照APG IV的科概念去整理family列。第三类是格式问题。全角空格、行尾多余的空格、半角括号混用、学名首字母没有大写每一样都会导致名称匹配失败。学名是拉丁语书写规范很严格属名首字母大写种加词全部小写中间用英文半角空格分隔。我在清洗数据时通常会写一小段R代码把字符串里的前后空格、多余空格全部去掉再检查一下是否有非ASCII字符混入spList$species - trimws(spList$species) spList$species - gsub(\\s, , spList$species) spList$genus - trimws(spList$genus) spList$family - trimws(spList$family)这三类坑补完之后你的名录才算真正达到可以交给V.PhyloMaker的状态。别嫌这一步麻烦实际经验是学名清洗花一小时后面能省出好几天。3. 核心实战从名录到系统发育树3.1 准备一个可以跑的示例数据为了让流程更直观我准备一个10个物种的示例数据。实际项目中你只需要把读入的名录整理成同样格式即可。spList - data.frame( species c( Fagus longipetiolata, Castanopsis eyrei, Quercus glauca, Machilus thunbergii, Cinnamomum camphora, Liquidambar formosana, Acer palmatum, Liriodendron chinense, Ulmus pumila, Ginkgo biloba ), genus c( Fagus, Castanopsis, Quercus, Machilus, Cinnamomum, Liquidambar, Acer, Liriodendron, Ulmus, Ginkgo ), family c( Fagaceae, Fagaceae, Fagaceae, Lauraceae, Lauraceae, Altingiaceae, Sapindaceae, Magnoliaceae, Ulmaceae, Ginkgoaceae ), stringsAsFactors FALSE )注意stringsAsFactors FALSE否则列会变成因子类型后面处理字符串时会遇到各种奇怪问题。如果你是从Excel读进来的数据读入之后最好先跑一遍str(spList)确认每列都是字符类型。3.2 调用phylo.maker核心函数就一个phylo.maker()。最基本的调用方式如下res - phylo.maker(spList, scenarios c(1, 2, 3))scenarios参数一次可以传多个值我习惯直接传c(1, 2, 3)这样能把三种不同处理策略的树一次性都生成出来后面做敏感性分析会很方便。运行结束后res是一个列表通常会包含tree、species.list、species.list.relative等几项。取树的时候我建议先看一下返回结构免得有的版本返回的是phylo对象、有的返回的是multiPhylo列表str(res) trees - res$tree if (inherits(trees, phylo)) { trees - list(trees) } tree_scen1 - trees[[1]] tree_scen2 - trees[[2]] tree_scen3 - trees[[3]]res$species.list这个表值得花时间认真看一遍。它记录了每个物种最终是怎么进入树的是本来就在超树中还是通过属插入还是通过科插入。这一张表就是你的匹配质量报告。如果通过科插入的物种太多说明名录里冷门种占比高后续结论需要谨慎解读。运行时间一般在几分钟以内主要取决于物种数量和名录的混乱程度。第一次加载超树会比较慢属正常现象。3.3 scenarios参数怎么选V.PhyloMaker的三个scenario是处理超树里查无此物的物种时的三种不同插接策略细微的拓扑差异集中在缺失物种的嫁接位置上。官方文档里有非常严格的算法描述这里我不逐条背书只说我实际使用时的判断逻辑如果你不确定该选哪一个就三个都跑一遍然后做敏感性分析。具体做法是分别用三棵树计算你要用的下游指标比如Faiths PD、MPD、MNTD然后看看结论是否一致。如果核心结论在三棵树之间是稳定的说明你对缺失物种的处理方式不敏感文章里可以直接报告其中一棵树的结果再在补充材料里说明三棵树结论一致。如果不同树之间指标差异很大说明名录里匹配失败的物种占比太高树本身的不确定性已经影响到研究结论了。这时候最优解不是纠结选哪个scenario而是回头去整理学名把匹配率提上去。我个人在常规项目中默认以scenario1那棵树为主分析树原因很简单它的插接策略在拓扑上相对保守不会对缺失物种在属内的精确位置做过多的假设推导。当然这只是个人使用习惯具体还需要结合你的科学问题来定。3.4 画出来看看树生成之后第一件事永远是可视化检查而不是直接丢进下游计算。我常用ggtree来画library(ggtree) library(ggplot2) p - ggtree(tree_scen1, layout fan, size 0.4) geom_tiplab(size 2.5, offset 0.2) xlim(-2, NA) p如果只关心拓扑结构和亲缘关系分组layout rectangular会更清晰。看到树的第一眼应该去确认几个关键节点你研究类群的主要支系是否聚在一起、科和属之间有没有出现明显异常的距离、是否有物种被孤零零挂在某个奇怪的位置。这个习惯能帮你发现spList里被忽略的分类学错误。想要进一步检查三棵树的差异可以用ape包做比较library(ape) all.equal.phylo(tree_scen1, tree_scen2, use.edge.length TRUE)差异通常只集中在少数嫁接物种身上如果差异面太大就要警惕是不是spList里存在系统性的学名错误。4. 拿到树之后质量检查与下游计算4.1 先检查树的几个基本属性拿到树就直接跑指标是我见过最多的翻车操作。至少要先做这三项基础体检library(ape) is.ultrametric(tree_scen1) is.binary(tree_scen1) Ntip(tree_scen1) Nnode(tree_scen1)is.ultrametric检查树是否是超度量树也就是所有末端到根节点的距离一致。很多生态学指标比如PD基于这一假设如果返回FALSE需要先处理。is.binary检查树是否是完全二叉的V.PhyloMaker在嫁接过程中有可能会产生软多分歧也就是某个节点下直接分出三个以上分支。Ntip应该等于你spList里的物种数如果不等说明有物种在匹配过程中被丢掉了要去检查res$species.list。4.2 处理多分歧如果is.binary(tree_scen1)返回FALSE最简单的做法是用ape::multi2di把多分歧随机解析成二叉set.seed(123) tree_binary - multi2di(tree_scen1, random TRUE)注意multi2di的随机解析会引入额外的拓扑不确定性。同一个多分歧不同的随机种子可能得到不同的二叉树而这会影响依赖最近邻居的一类指标比如MNTD。建议的做法是设置固定随机种子保证结果可复现更严谨一点可以生成多棵随机解析树分别计算指标后取平均值和范围。前者适合日常探索后者适合正式发表。4.3 对接picante计算常见指标树和样本数据都准备好之后就可以接上生态学里常用的R包picante了。这里最常见的一个坑是样地-物种矩阵的列名必须和树的tip.label完全一致包括顺序、大小写、空格都不能差。我每次都先用match.phylo.comm对齐一下library(picante) comm - matrix( c( 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1 ), nrow 2, byrow TRUE, dimnames list( c(Plot_A, Plot_B), spList$species ) ) match_result - match.phylo.comm(tree_scen1, comm) phylo_comm - match_result$comm对齐之后就可以放心计算了pd_out - pd(phylo_comm, tree_scen1, include.root TRUE) mpd_out - mpd(phylo_comm, cophenetic(tree_scen1)) mntd_out - mntd(phylo_comm, cophenetic(tree_scen1))pd得到每个样地的Faiths PD和物种丰富度mpd和mntd分别得到平均成对距离和平均最近邻距离。这几个指标是群落系统发育分析里的老面孔了。再往后的标准化效应大小比如NRI和NTI可以基于空模型做随机化得到逻辑也是在这个基础之上。5. 常见问题与避坑实录5.1 匹配率太低怎么办匹配率低的头号原因是学名不一致。我处理过的最典型情况名录里写的是Quercus glauca但实际研究区域内这个种在最新分类系统中被划到了别的属或者名录里用的是中国植物志的老科名而超树用的是APG IV体系。遇到这种情况先不要急着怪工具回到名录本身去校对。校对这一步我强烈建议在Excel里做一次人工抽查随机抽20个物种去POWO上核对有效名如果20个里有3个以上对不上说明整个名录都要系统地清理一遍。5.2 科名过时或不一致超树执行名称匹配的时候科名是最后一道兜底防线。如果科名对不上这个物种就会变成彻底无法定位的状态。注意一个容易忽略的细节有些属名在历史上曾经被放在不同的科下同样是Liquidambar formosana老资料里写Hamamelidaceae新系统写Altingiaceae如果你名录里混用了两套科名即便物种名本身正确也可能在匹配过程中出问题。所以我在准备spList之前会专门写一个函数把科名校验到APG IV框架下统一后再交给phylo.maker。5.3 后续分析报错的几个典型场景很多人反应pd()或mpd()报错仔细一看都是树的问题有的是树里有多分歧导致cophenetic生成了Inf有的是样地矩阵的列名和tip.label不匹配还有的是树不是ultrametricpd直接给出NA。我的建议是每一步都做一个轻量确认算之前先all(sort(colnames(comm)) sort(tree$tip.label))再is.ultrametric(tree)都通过了再进分析。这套检查看起来啰嗦但能挡住大多数莫名其妙的报错。5.4 常见问题速查表症状可能原因解决建议大量物种未匹配名录含异名、拼写错误或不规范命名用POWO等权威名录校对有效名提示找不到科科名使用了老分类系统统一到APG IV科概念运行时报内存不足内置超树较大内存消耗高重启R关闭其他程序至少8GB内存输出树含有软多分歧缺失物种嫁接导致用multi2di随机解析设置随机种子pd或mpd报错树非ultrametric或含多分歧先is.ultrametric、is.binary检查comm和树匹配报错样地矩阵列名与tip.label不一致用match.phylo.comm对齐不同scenario结果差异明显名录中匹配失败的物种占比较高优先做学名校正并报告敏感性分析结果我在实际项目里有一个固定的处理策略先跑phylo.maker(spList, scenarios c(1, 2, 3))然后看res$species.list里的匹配统计。如果通过科插入的比例超过10%我会再花一天时间回去校对学名。校对完再跑一遍通常匹配率能明显改善。如果重新梳理之后仍然有比较高的插入比例那就说明这个类群的超树本身覆盖不足必须在论文里如实报告这一限制。这个操作习惯不算复杂但能避免很多后续审稿人提出你的树可信度如何时的被动局面。