用PartitionFinder2 筛选系统发育分析分隔模型

PartitionFinder2是用来在建立进化树时筛选最优进化模型以及分区（分隔）方案的软件，是用Python2语言编写的。

核苷酸序列在不同位点有不同的突变速率。核苷酸序列又分为编码基因和非编码基因。编码基因中，密码子第一、二位往往较为稳定，第三位往往变异速率较高。非编码基因因为受到的选择压力一般较小，所以往往可保留更多突变。不同基因以及不同位点的突变速率不同，可能对所推断进化树的稳定性有很大影响。所以在多基因建立进化树过程中，设置分隔模型就显得很重要。

但是分隔模型怎样设置才算合理呢？在PartitionFinder软件开发出来之前，研究人员一般通过将不同基因分开，例如gene1、gene2、gene3，再将每个基因的不同位点分开如gene1_1、gene1_2、gene1_3作为分隔模型，以提高进化树的稳健性。

RAxML、MrBayes、BEAST等常用系统发育软件都支持分隔模型，但并不能用于确定最优化的分隔模型设定方案。在建树过程中，如果设置的分隔模型过多，则拟合的参数会过多，造成结果不准确；如果设置的分隔模型过少，或者不设置分隔模型、设置的不合理，也会造成进化树不准确。很多学者已经意识到这个问题，但是一直苦于没有很好的应对方法。

赤池信息量AIC整合了模型的似然值（Likelihood）以及所要估计的参数数量。模型选择的简约理论认为，AIC或者BIC最小的模型是最优的模型。如果似然值已经足够高，进一步增加模型的参数已经不能再显著提高似然值的情况下，可以认为已经找到了最优模型。但是，对于分隔模型的筛选来说，碱基替换模型与分区方式的组合数量非常多，如果要对每一种分区方案对应的碱基模型筛选出最优组合，在很多情况下几乎是不可能的。这是因为：首先，从碱基比对矩阵计算Likelihood是十分耗费时间的。其次，分隔模型的各种组合的数量呈几何级数增长，如果通过几个基因建树，则可能的分隔模型的数量已经超出了大部分计算机的计算能力。此时就需要引入分隔模型的启发式搜索（Heuristic Search）。

PartitionFinder的作者从原理上解决了以上的问题，并通过Python语言实现了相应的算法（Lanfear et al., 2012）。作者证明，PartitionFinder所得的分隔模型比之前的简单处理更加合理。不仅如此，PartitionFinder在获得最优化分隔方案后，同时会给出每个分隔模块所对应的最优进化模型。因此，Modeltest、jModeltest以及ProtTest等软件都已经被PartitionFinder所超越。事实上，2012年Lanfear介绍PartitionFinder的论文发表后，在google scholar上已经被引用了3223次（截至2019年1月19日）。

本文简述PartitionFinder的安装和使用，供同行参考。

1. 下载和安装

1.1 下载

PartitionFinder的下载网址http://www.robertlanfear.com/partitionfinder/，下载后用WinRAR或者7zip解压缩。

图1 PartitionFinder的结构

其中:

  • docs为说明文档
  • Examples 文件夹下为示例文件，包含DNA序列、形态以及氨基酸序列三个例子
  • partfinder 文件夹下为实现算法的 python脚本
  • Programs 文件夹下为PHYML可执行文件，用来计算likelihood
  • Submodules 文件夹下包含了raxml文件夹，但是该文件夹为空。作者并未交代这个文件夹的内容。
  • tests 文件夹下为开发时的测试文件
  • PartitionFinder.py 为筛选DNA序列最优分隔模型的Python脚本
  • PartitionFinderMorphology.py 为筛选形态数据最优分隔模型的Python脚本
  • PartitionFinderProtein.py 为筛选氨基酸序列最优分隔模型的Python脚本

运行PartitionFinder时，要根据数据类型分别调用三个Python脚本，不能混淆。

1.2. 安装Anaconda

PartitionFinder是用Python2开发的，依赖于numpy、pandas、pytables、pyparsing、scipy、sklearn 等若干程序包。在Windows下，安装这些程序包往往较为麻烦。作者Lanfear推荐已经预装了这些程序包的Anaconda Python。

Anaconda Python的下载地址为：

https://www.anaconda.com/distribution/

PartitionFinder只能在Python2.7上运行，不能使用Python3。因此，这里只能下载Anaconda Python2.7。

双击开始安装。

选择将以下路径添加到系统路径:

C:\ProgramData\Anaconda2;C:\ProgramData\Anaconda2\Scripts; C:\ProgramData\Anaconda2\Library\bin; 以便程序能从命令行CMD运行。

图2 下载Anaconda Python2.7

2. 所需文件

运行PartitionFinder，用户需要提供：

  1. phylip文件：包含所要读取的序列
  2. partition_finder.cfg文件：包含模型筛选时设定的参数

2.1 PHYLIP文件

phylip文件可以通过MUSCLE、CLUSTAL、MAFFT等比对序列的软件生成。如果是多基因序列，一般对每个基因分别比对，手工调整之后，再通过Geneious或phylotools（https://github.com/helixcn/phylotools）等软件生成supermatrix，以phylip格式导出即可。

PHYLIP格式的详细定义，参考

http://www.atgc-montpellier.fr/phyml/usersguide.php?type=phylip。

PartitionFinder可以接受relaxed phylip，文件中的物种名最长可以达100个字符。

图3 Phylip文件

2.2 partition_finder.cfg 文件

PartitionFinder可以处理三种类型的数据，分别为：

  1. 编码区DNA序列
  2. 氨基酸序列或非编码区DNA序列
  3. 形态数据

根据数据类型的差别，partition_finder.cfg文件会各有不同。

  1. 编码区DNA序列的配置文件要提供每个片段的
  3. 起始和终止位置;
  4. 编码位置（1、2、3）
  5. 氨基酸序列和非编码区DNA序列，只需提供每个基因片段的起始和终止位置。
  6. 形态数据，无需直接指定模块，在计算时采用k-means聚类等算法以划分成不同的组别。

图4 partition_finder.cfg文件的内容

3. 运行

对于DNA序列：

  1. 在partitionFinder文件夹下创建一个新文件夹，命名为test20190219

  2. 将examples/nucleotide文件夹下的test.phy文件以及partition_finder.cfg配置文件拷贝到test20190219文件夹下

  3. 在CMD中，将工作目录转换到partitionFinder下

  4. 调用PartitionFinder.py，格式如下：

python PartitionFinder.py test20190219

  1. 生成的结果会自己保存在test20190219文件夹下，并生成若干新文件（如log.txt）以及文件夹analysis。

  2. 模型筛选的结果就保存在analysis/best_scheme.txt文件中。

图5 运行partitionFinder

图6 partitionFinder的运行结果

4. cfg配置文件

数据类型不同的PHYLIP文件，需要在partition_finder.cfg文件中设置不同的参数。

4.1 cfg文件的参数

4.1.1 alignment

指出phylip的文件名

1
alignment = test.phy;

4.1.2 branchlengths

在优化进化树枝长用于Likelihood计算时，是否各枝长一起优化。

1
branchlengths = linked;

4.1.3 models

models指定模型筛选的范围

1
models = all;

4.1.4 model_selection

model_selection模型筛选的标准，一般用aicc

1
model_selection = aicc;

4.1.5 data_blocks

datablocks 设定数据分区，例如：

1
2
3
4
[data_blocks]
Gene1_pos1 = 1-789\3;
Gene1_pos2 = 2-789\3;
Gene1_pos3 = 3-789\3;

模型筛选的算法 search = greedy;

4.2 部分参数详解

4.2.1 branchlengths

枝长是否相互关联，即是否对每个分隔模式独立估计进化树的枝长（目前只有BEAST、MrBayes、RAxML支持独立估算枝长）

4.2.2 不同数据类型模型筛选的范围

默认情况下models=all；如果用户需要使用更多模型，可设定models=allx，但应尽量避免，这是因为对少数几个模型进行优化要耗费很多时间。

  1. 对于离散性状，例如01性状（BINARY）以及无序质量性状，模型为(MULTISTATE)BINARY+G、BINARY+G+A、MULTISTATE+G、MULTISTATE+G+A，相关详情可参考RAxML的说明书。注意：性状数据模型不能计算AIC。
  2. 对于DNA****序列，models=all将比较内置的56种模型。即由14种基本模型为基础：JC、K80、TrNef、K81、TVMef、TIMef、SYM、F81、HKY、TrN、K81uf、TVM、TIM、GTR，同时参考是否用Gamma或者I来表征不同位点进化速率的差异。注意：在分析DNA数据时，如果添加了 --raxml选项，则只计算GTR、GTR+G、GTR+I+G三种模型的AIC的值。
  3. 对于氨基酸序列，默认情况下（使用PHYML计算），要比较的基本模型为：LG、WAG、MTREV、DAYHOFF、DCMUT、JTT、VT、BLOSUM62、CPREV、RTREV、MTMAM、MTART、HIVB、HIVW，每种又分成8种情况：LG、LG+F、LG+G、LG+G+F、LG+I、LG+I+F、LG+I+G、LG+I+G+F ，其中F为是否使用经验氨基酸频率。如果开启了RAxML选项，则会比较LG、WAG、MTREV、DAYHOFF、DCMUT、JTT、VT、BLOSUM62、CPREV、RTREV、MTMAM、MTART、HIVB、HIVW、MTZOA、PMB、JTTDCMUT、FLU、STMTREV、DUMMY、DUMMY2等21种模型，每个模型考虑6种情况。

所筛选模型范围，如果改为models = allx，将进一步扩大：

  1. 对于DNA序列，使用PHYML时，会比较84种模型；使用RAxML时，会比较6种模型；
  2. 对于氨基酸序列，使用PHYML不能获得结果，所以只能使用RAxML，比较的模型数量为195个。此时，所有的碱基频率或者氨基酸频率都是用极大似然估计获得的。

若从指定的模型中筛选：

  1. 如果后续用mrbayes或者beast做分析，则需要设定models为 models = mrbayes;或 models = beast;，以便只从这两个软件中有的模型中筛选。
  2. 用户也可设定PartitionFinder只考虑有Gamma或Gamma+I的相关的模型
1
2
models = gamma;
models = gammai;
  1. 或直接指定DNA模型的筛选范围:
1
models = JC, JC+G, HKY, HKY+G, GTR, GTR+G;
  1. 指定氨基酸模型的筛选范围
1
models = LG, LG+G, LG+G+F, WAG, WAG+G, WAG+G+F;
  1. 指定模型选择的标准，作者推荐使用corrected AIC, 即AICc
1
2
# model_selection: AIC | AICc | BIC
model_selection = aic;

4.2.3 设定分区

分区模块需要用[data_blocks]标记。

1
2
3
4
5
[data_blocks]
Gene1_codon1 = 1-1000\3; 
Gene1_codon2 = 2-1000\3;
Gene1_codon3 = 3-1000\3;
intron = 1001-2000;

[data_blocks] 后紧跟基因以及位点的名称。如上面代码所示： 数据包括Gene1以及intron两部分。对于编码基因Gene1，由于密码子的1、2、3位的进化速率可能不同，因此可将1-1000位点的1、2、3碱基分开，分别作为单独的分区。

上述命令中：

Gene1_codon1 = 1-1000\3; 表示：在第1-1000位点，从第1个碱基开始，步长为3的所有碱基

Gene1_codon2 = 2-1000\3; 表示：在第2-1000位点，从第2个碱基开始，步长为3的所有碱基

Gene1_codon3 = 3-1000\3; 表示：在第3-1000位点，从第3个碱基开始，步长为3的所有碱基

intron = 1001-2000; 表示：从第1001-2000位点。由于是内含子，没有编码功能，所以不用将1、2、3位碱基分开。

如果是为MrBayes筛选分区模型，可在每行前增加charset关键词。

1
2
3
4
charset Gene1_codon1 = 1-1000\3;
charset Gene1_codon2 = 2-1000\3;
charset Gene1_codon3 = 3-1000\3;
charset intron = 1001-2000;

4.2.4 搜寻策略

根据数据量以及要筛选的分区组合的数目，PartitionFinder的搜索策略，要通过修改search参数做相应调整。

1
search: all | greedy | rcluster | rclusterf | hcluster | kmeans | user

其中:

  • search = all 一般用于12个分区以下；
  • search = greedy 10-100个分区；
  • search = rcluster 100个以上的分区；
  • search = rclusterf 只适用于RAxML选项开启时，尤其适用于有若干分析需要数据很大，同时ML需要优化很长时间时，可减少等待的时间；
  • search = hcluster 在研究中不推荐使用，仅为了开发而保留；
  • search = user 如果用户要对上述定义的分区组合，按照指定的方式筛选， 则可在search = user命令之后，指定搜索的策略。同一对括号中的基因表示将作为一个整体进行分析（参见英文版说明书第21页）
1
2
3
4
together= (Gene1_codon1, Gene1_codon2, Gene1_codon3, intron); 
intron_123  = (Gene1_codon1, Gene1_codon2, Gene1_codon3) (intron); 
intron_12_3 = (Gene1_codon1, Gene1_codon2) (Gene1_codon3) (intron); 
separate= (Gene1_codon1) (Gene1_codon2) (Gene1_codon3) (intron);

4.2.5 搜寻策略调整举例

4.2.5.1 少于12个分区

1
2
3
4
branchlengths = linked;
models = all;
model_selection = aicc;
search = all;

4.2.5.2 对10-100个分区进行筛选

超过100个等位基因的序列，需要用贪婪算法，基本步骤如下：

  1. 设定每个基因的起始位点以及结束位点，编码位置（1，2，3）等
  2. 设定.cfg文件
  3. 用raxml计算Likelihood

对于DNA为：

python PartitionFinder.py InputFoldername --raxml

对于氨基酸为：

python PartitionFinderProtein.py InputFoldername --raxml

4.2.5.3 超大数据（1000个）的分区筛选

超大型数据的位点（1000个以上），贪婪算法仍然太慢，此时推荐使用其他算法，以筛选最优分隔模型设置方案（Lanfear et al 2014），基本步骤如下:

  1. 设定每个基因的起始位点以及结束位点，编码位置（1、2、3）等
  2. 更改.cfg文件，设定 search = rcluster
  3. 在cmd中运行PartitionFinder:
1
python PartitionFinder.py InputFoldername --raxml

如仍然过慢，则更改recluster-max 参数:

1
python PartitionFinder.py InputFoldername --raxml --rcluster-max 100

4.2.6 用户自定义进化树

在模拟研究中，可能会假设已知进化树的结构，探讨能否筛选到最优分区模型。此时可指定一棵进化树，在后续的分析中，该进化树的拓扑结构将被保留，但是枝长会根据GTR+I+GAMMA模型重新估计。指定进化树的命令如下：

1
2
3
ALIGNMENT FILE # 
alignment = test.phy; 
user_tree_topology = tree.phy

5.输出结果

生成的结果位于同一文件夹下的analysis文件夹中，该文件夹包含以下文件：

  1. best_schemes.txt纯文本文件，包括所筛选到的最优化分隔模型划分方案以及每个区块对应的最优模型以及RAxML和Nexus的分隔模型。
  2. subsets文件夹:每个分区所对应的AICc值以及相应的筛选过程
  3. schemes文件夹:所分析的分隔模型划分方案，详细的分析过程

6.命令行选项

在运行PartitionFinder时，可在命令行设定一些参数，具体如下：

  • --all-states 仅用于k-means选项开启的情况，用于限定状态的数量
  • --force-restart 删除之前的分析，完全重新开始
  • --min-subset-size 设定子集的大小。只用于k-means算法（这个选项在说明书中交代得不是很清楚）
  • --no-ml-tree 设定此选项后，PartitionFinder将从Neighbour Joining（PHYML）或者Parsimony树（RAxML）开始。
  • --processors N, -p N 多线程计算，使用多个CPU内核。默认情况下，PartitionFinder会使用多个核，如果要控制所用到的核数，则应做相应设定。
  • --quick, -q 让PartitionFinder停止写一些记录文件，如果处理的数据特别大，有可能会提升速度，常用于使用贪婪算法的情况下。
  • --raxml 使用RAxML帮忙寻找最优分隔模型。
  • --rcluster-max N 设定cluster的数量，默认值N为1000或者分区数乘以10中的较大者，一般已经可以保证搜寻到正确的分隔模型。
  • --rcluster-percent N 设定要搜寻的所有可能的分区组成的百分数，默认为10%。该参数配合--rcluster-max N，可减少重复计算。
  • --save-phylofiles 保存数据分析过程中生成的进化树，主要用来纠错。
  • --weights "W_rate , W_base , W_model , W_alpha" 只用于--raxml或hcluster or rcluster选项开启时，分别表示几个参数的权重（参考说明书第26页）:
1
2
3
4
1. the overall rate for a subset
2. the base/amino acid frequencies
3. the model parameters
4. the alpha parameter (which describes  gamma distributed rates across sites)

进一步阅读

  • Lanfear, R., Calcott, B., Ho, S. Y., & Guindon, S. (2012). PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Molecular biology and evolution, 29(6), 1695-1701.

附录

DNA序列碱基替换模型的三大要素

1. 碱基频率的估算

模型需要提供ATCG频率的估计方法，目前有四种方法:

  1. equal 直接假定ATCG的频率相同
  2. model 是从其他数据读取不同氨基酸的频率（仅用于氨基酸序列）
  3. empirical 是直接从数据计算ATCG的频率
  4. ML是用极大似然法估计ATCG的频率

注意：在模型后添加+X，即可开启ML方法估计ATCG频率，例如F81+X，不过一般用户都是用empirical的方法估计频率，而且结果与ML的结果相差无几。

2. 碱基替换速率矩阵

用来表征碱基之间的转换速率（Relative rates of substitution），其中GTR模型最为复杂，需要估计所有6个参数。

3. 不同位点进化速率的差异

分成以下几种情形：

  1. 各位点速率恒定；
  2. 部分位点进化树进化树速率恒定不变 （I用来表示进化速率恒定位点的比例）：GTR + I；
  3. 假设各位点的进化速率符合Gamma分布（因为Gamma分布要估计的参数少，同时曲线有足够的变化，能够较好地拟合多数情况，所以用Gamma分布）：GTR + GAMMA；
  4. 同时考虑各位点的进化树速率的差异以及恒定位点的比例：GTR + GAMMA + I。