一步一步教你构建进化树（ML树）
大多数人习惯了利用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进行展示
欢迎关注生信人