大多数人习惯了利用MEGA构建NJ树，速度快，准确性也不是很差。但是比较严格做法是，构建多种树进行比较，比较常用的那就是ML树。高质量期刊一般会采用NJ树和ML树相互验证的方法。构建ML树，速度较慢，大家需要注意。构树常用的软件流程基本上可以采用MUSCLE+PhyML了。

1. MUSCLE

由于利用PhyML构建ML树，需要phy格式的比对文件，因此需要用MUSCLE产生。需要注意的一点是phy格式比对文件对序列ID要求最多10个字符，因此构树之前要进行更改。如下图，Glyma.10G0这个ID如果不更改，最终比对完成后会显示ID的前十个字符，其全部的ID是Glyma.10G010000。

修改的脚本参照下面：

my $num=1;

`mkdir $od` unless (-d "$od");

$od=“Change_ID”###输出目录，可以自定义

 $fa=“test.fa”;###序列文件名字，改成你自己的即可

$index="mapk";###4个字符以内的任意前缀

my %VS;

open (OUT,">$od/VS.txt")|| die "cannot open $od:$!";###ID 改前后对照表

open (SEQ,">$od/change_name.fa")|| die "cannot open $od/change_name.fa:$!";

my $ina = Bio::SeqIO->new(-file => $fa, -format => 'fasta');

while(my $obj = $ina->next_seq()){

        my $id = $obj->id;

        my $seq = $obj->seq;

        my  $id2="$index$num";

         $VS{$id}=$id2;

         $num++;

         print SEQ ">$id2\n$seq\n";

}

close SEQ;

foreach my $key(keys %VS){

        print OUT "$key\t$VS{$key}\n";

}

修改完成后进行MUSCLE比对，注意比对输出格式选择：

muscle -in change_name.fa -phyiout change_name.fa.phy

2. ML树构建。

因为ML构建速度非常慢，所以PhyML支持多线程，可以通过设置多个CPU解决速度慢的问题（需要openmpi软件支持）。使用的示例命令如下：

openmpi/bin/mpirun -n  10  phyml/20151210/bin/phyml-mpi -i change_name.fa.phy  -d nt -b 1000

其中主要参数意思如下：

-n：CPU数目

-i: 比对文件

-d：序列类型，核酸选择nt，蛋白选择aa

-b：bootstrap次数，一般设置1000

最终生成的树文件为 change_name.fa.phy_phyml_tree。

3. Bootstrap值与ID转换。

生成的nwk格式的进化树。

(mapk14:0.14773305,(mapk15:0.08940692,mapk5:0.09227976)1000:0.16997388,(mapk18:0.37695718,(mapk6:0.62762808,(((mapk9:0.02521634,mapk13:0.02897656)1000:0.09387741,(mapk8:0.06673311,mapk3:0.06951497)702:0.05152950)1000:0.56517119,((mapk10:0.10842097,(mapk17:0.10174175,(mapk7:0.06451028,mapk2:0.08751185)905:0.06671077)650:0.01753958)630:0.04626638,((mapk16:0.06591709,(mapk4:0.05730643,mapk1:0.07344786)909:0.05745075)890:0.00935337,(mapk12:0.02930148,mapk11:0.02615966)1000:0.06766173)701:0.05070130)600:0.13183564)1000:0.72763617)940:0.19044937)1000:0.57474285);

可以看到bootstrap值为具体的值，而不是一个百分比例（举个例子，1000次有650次能够支持，比例应为65，所以需要将650改成65）。

下面脚本进行Bootstrap和ID转换。

$out="final.tree";

$index="mapk";

open(OUT,">$out");

my %VS1;

open (IN,"$VS")|| die "cannot open $VS:$!";##$VS改成步产生的VS文件，改回原先的ID

while(<IN>){

chomp;

my @tmp=split/\s+/;

$VS1{$tmp[1]}=$tmp[0];

}

my %VS2;

open (TREE,"$tree")|| die "cannot open :$!";###$tree改成生成的tree文件

while(<TREE>){

my @tmp=split/\:/;

for(my $i=0;$i<=$#tmp;$i++){

        if($tmp[$i]=~/^(.*\))(\d+)$/){

                my $ratio=int($2/$num*100);

                $tmp[$i]="$1$ratio";

        }

}

for(my $i=0;$i<=$#tmp;$i++){

        if($tmp[$i]=~/^(.*)($index\d+)$/){

                my $name=$VS1{$2};

                $tmp[$i]="$1$name";

        }

}

print OUT join(":",@tmp)."\n";

}

4.树展示

利用MEGA中的User Tree中的Display Newick Trees进行展示

欢迎关注