使用 PhyML 构建进化树

1. PhyML 简介

使用 PhyML 构建最大似然树。
参考文献:New Algorithms and Methods to Estimate Maximum-Likelihood Phylogenies: Assessing the Performance of PhyML 3.0

2. PhyML 的下载和安装

$ wget http://www.atgc-montpellier.fr/download/binaries/phyml/PhyML-3.1.zip
$ unzip PhyML-3.1.zip
$ mv PhyML-3.1 /opt/biosoft/
$ ln -s /opt/biosoft/PhyML-3.1/PhyML-3.1_linux64 /opt/biosoft/PhyML-3.1/PhyML
$ echo 'PATH=$PATH:/opt/biosoft/PhyML-3.1/' >> ~/.
$ source ~/.bashrc

3. PhyML 的使用

PhyML 的输入文件为 phylip 格式。
常用例子:

$ PhyML -i proteins.phy -d aa -b 1000 -m LG -f m -v e -a e -o tlr

常用参数:

-i seq_file_name
输入文件,phylip 格式的多序列比对结果。
-d data_type default:nt
该参数的值为 nt, aa 或 generic。
-b int
设置 bootstrap 次数。
-m model
设置替代模型。 核酸的模型有: HKY85(默认的), JC69, K80, F81, TN93, GTR ; 氨基酸的模型有:LG (默认的), WAG, JTT, MtREV, Dayhoff, DCMut, RtREV, CpREV, VT, Blosum62, MtMam, HIVw, HIVb 。
-f e,m or fA,fC,fG,fT
设置频率计算的方法。 e 表示使用比对结果中不同氨基酸或碱基出现的频率来计算; m 表示使用最大似然法计算碱基频率,或使用替换模型计算氨基酸频率; fA,fC,fG,fT 则是 4 个浮点数,表示 4 中碱基的频率,仅适合核酸序列。
-v prop_invar
设置不变位点的比例,是一个[0,1]区间的值。或者使用 e 表示程序获得其最大似然估计值。
-a gamma
gamma 分布的参数。此参数值是个正数,或者使用 e 表示程序获得其最大似然估计值。在 ProtTest 软件给出的最优模型中含有 G 时,使用该参数。
-o params
参数优化的选项。t 表示对 tree topology 进行优化; l 表示对 branch length 进行优化; r 表示对 rate parameters 优化。
params=tlr 这表示对 3 者都进行优化。 params=n 表示不进行优化。

4. PhyML 结果

PhyML 的输出结果为:

proteins.phy_phyml_tree.txt        :    最大似然法构建的进化树
proteins.phy_phyml_boot_stats.txt  :    bootstrap 的统计信息
proteins.phy_phyml_boot_trees.txt  :    bootstrap 树
proteins.phy_phyml_stats.txt       :    程序运行的中的参数和结果统计

使用 ProtTest 来选择最优氨基酸替代模型

1. ProtTest 简介

ProtTest 用来进行最优氨基酸替代模型的选择。相应的,适用于核苷酸的软件是 jModeltest。
ProtTest 通过 PhyML 对进化树和模型参数的最大似然估计,通过 AIC, BIC 分值或 DT 来寻找最佳模型。分值越小越优。
ProTest 3.2 版本包含 15 种不同类型的 rate matrices;考虑到位点的 rate variation (+I: invariable sites; +G: gamma-distributed rates) 和 observed amino acid frequencies (+F), 共有 120 种不同的模型。
ProtTest 官网:https://code.google.com/p/prottest3/
从此处下载该软件。可能需要设置代理后下载。
参考文献:ProtTest 3: fast selection of best-fit models of protein evolution

2. ProtTest 下载和安装

$ tar zxf prottest-3.4-20140123.tar.gz -C /opt/biosoft
$ cd /opt/biosoft/prottest-3.4-20140123
$ echo 'export PROTTEST_HOME=/opt/biosoft/prottest-3.4-20140123' >> ~/.bashrc
查看说明文档:
$ less README

3. ProtTest 的使用

ProtTest 使用 JAVA 编写,有图形化和命令行两种运行模式。

3.1 图形化界面使用

必须要进入到程序的所在的目录运行程序以启动图形化界面
$ cd /opt/biosoft/prottest-3.4-20140123/runXProtTestHPC.sh
$ runXProtTestHPC.sh

启动 JAVA 界面后,点击 File–Load Alignment, 上传多序列比对结果;然后点击 Analysis–Compute likehood scores, 选择所使用的线程数,以及候选模型的选择,和计算 likelihood 的 topology;然后点击 Compute, 进行计算,所需要消耗的实际有点长;计算完毕后,点击 Selection–Results 来查看结果。通过 AIC, BIC, AICc 和 DT 来查看其得分,点击表格的第1行进行排序,寻找分值最小的模型作为最优氨基酸替代模型。

3.2 命令行运行

常用例子:

不加参数运行,则给出帮助文档:
java -jar /opt/biosoft/prottest-3.4-20140123/prottest-3.4.jar

常用的命令行:
java -jar /opt/biosoft/prottest-3.4-20140123/prottest-3.4.jar -i proteins.phy -all-distributions -F -AIC -BIC -tc 0.5 -threads 24 -o prottest.out

ProtTest 的常用参数:

-i alignment_filename
必须参数,输入多序列比对结果文件。
-o output_filename
输出的文件名。不设置,则默认输出到标准输出。
-[matrix]
指定需要分析的 matrix 。 该 matrix 可以被替换为 JTT LG DCMut MtREV MtMam MtArt Dayhoff WAG RtREV CpREV Blosum62 VT HIVb HIVw FLU 这 15 种 matrix。 若不指定,则默认全选。
-all-distributions
指定 matrix 模型结合 I 或 G 或 I+G
-F
指定 matrix 模型结合 empirical frenquency estimation
-AIC
输出结果中按 AIC (Akaike Information Criterion) 排序
-BIC
输出结果中按 BIC (Bayesian Information Criterion) 排序
-AICC
输出结果中按 AICc (Corrected Akaike Information Criterion) 排序
-DT
输出结果中按 DT (Decision Theory Criterion) 排序
-tc consensus_threshold
输出满足指定阈值的一致树。该值在 0.5~ 1.0 之间。[0.5 = majority rule consensus ; 1.0 = strict consensus]
-threads number_of_threads
使用的 CPU 数。

使用 Gblocks 提取保守序列

1. Gblocks 简介

Gblocks用于从多序列比对结果中提取保守位点,以利于下一步的进化分析。
在线说明文档:http://molevol.cmima.csic.es/castresana/Gblocks/Gblocks_documentation.html
在线服务网站:http://molevol.cmima.csic.es/castresana/Gblocks_server.html

2. Gblocks 安装

$ wget http://molevol.cmima.csic.es/castresana/Gblocks/Gblocks_Linux64_0.91b.tar.Z
$ sudo yum install -y ncompress
$ tar Zxf Gblocks_Linux64_0.91b.tar.Z -C /opt/biosoft/
$ echo 'PATH=$PATH:/opt/biosoft/Gblocks_0.91b/' >> ~/.bashrc
$ source ~/.bashrc

3. Gblocks 使用

Gblosk 有两种使用方式,第一种是交互式的方式(按提示输入文件改变参数),第二种是命令行式(在命令行中输入参数)。
命令行式的常用例子:

使用默认的设置:
$ Gblocks proteins.fasta -t=p
必须是 fasta 文件在前,参数在后。若没有参数,则进入交互式界面。

得到更长的序列:
$ Gblocks proteins.fasta -b4=5 -b5=h

命令行式的常用参数:

-t= default:p
设置序列的类型,可选的值是 p,d,c 分别代表 protein, DNA, Codons 。
-b1= default:( 序列条数的 50% + 1 )
设定保守性位点必须有 >= 该值的序列数。该参数后接一个 integer 数,默认下比序列条数的 50% 大 1.
-b2= default: 序列条数的 85%
确定保守位点的侧翼位点时,其位点必须有 >= 该值的序列数。该值必须要比 -b1 的值要大。
-b3= default: 8
最大连续非保守位点的长度。
-b4= default: 10
保守位点区块的最小长度。该值必须 >=2 。
-b5= default: n
设置允许含有 Gap 位点。可选的值有 n,h,a 分别代表 None, With Half, All 。 当为 h 时,表示
-e= default: -gb
设置输出结果的后缀。

使用 MAFFT 进行多序列比对

1. MAFFT 简介

最经典和广为熟知的多序列比对软件是 clustalw 。 但是现有的多序列比对软件较多,有文献报道:比对速度(Muscle>MAFFT>ClustalW>T-Coffee),比对准确性(MAFFT>Muscle>T-Coffee>ClustalW)。因此,推荐使用 MAFFT 软件进行多序列比对。

2. MAFFT 下载安装

$ wget http://mafft.cbrc.jp/alignment/software/mafft-7.158-without-extensions-src.tgz
$ tar zxf mafft-7.158-without-extensions-src.tgz
$ cd mafft-7.158-without-extensions/core
$ perl -p -i -e 's#PREFIX =.*#PREFIX = /opt/biosoft/mafft#' Makefile
$ perl -p -i -e 's#BINDIR =.*#BINDIR = /opt/biosoft/mafft/bin/#' Makefile
$ make
$ make install
$ echo 'PATH=$PATH:/opt/biosoft/mafft/bin/' >> ~/.bashrc
$ source ~/.bashrc

检测软件是否正确安装
$ cd ../test
$ rehash                                                   # if necessary
$ mafft sample > test.fftns2                               # FFT-NS-2
$ mafft --maxiterate 100  sample > test.fftnsi             # FFT-NS-i
$ mafft --globalpair sample > test.gins1                   # G-INS-1
$ mafft --globalpair --maxiterate 100  sample > test.ginsi # G-INS-i
$ mafft --localpair sample > test.lins1                    # L-INS-1
$ mafft --localpair --maxiterate 100  sample > test.linsi  # L-INS-i
$ diff test.fftns2 sample.fftns2
$ diff test.fftnsi sample.fftnsi
$ diff test.gins1 sample.gins1
$ diff test.ginsi sample.ginsi
$ diff test.lins1 sample.lins1
若 diff 的结果不换回异常,则正确安装。

3. MAFFT 使用

MAFFT 有一些参数,适合不同情况下的多序列比对。
软件输入为 FASTA 格式,默认输出到标准输出。

3.1 较为精确的方法

L-INS-i

最准确的方法。适合于 <200 条序列,且序列长度 <~2000 aa/nt 的比对。

$ mafft –localpair –maxiterate 1000 input [> output]
$ linsi input [> output]

G-INS-i

适合于序列长度相似的多序列比对。序列条数 <200, 序列长度 <~2000 aa/nt 。

$ mafft –globalpair –maxiterate 1000 input [> output]
$ ginsi input [> output]

E-INS-i

适合序列中包含较大的非匹配区域。序列条数 <200, 序列长度 <~2000 aa/nt 。

$ mafft –ep 0 –genafpair –maxiterate 1000 input [> output]
$ einsi input [> output]

3.2 节约速度的方法

FFT-NS-i

减少迭代次数,最大迭代次数减为 2 。

$ mafft --retree 2 --maxiterate 2 input [> output]
$ fftnsi input [> output]

FFT-NS-2

最大迭代次数减为 0 。

$ mafft --retree 2 --maxiterate 0 input [> output]
$ fftns input [> output]

FFT-NS-1

此方法非常快速,适合 >2000 条序列的多序列比对。

$ mafft --retree 1 --maxiterate 0 input [> output]

NW-NS-i

迭代过程中不进行 FFT aproximation

$ mafft --retree 2 --maxiterate 2 --nofft input [> output]
$ nwnsi input [> output]

NW-NS-2

$ mafft --retree 2 --maxiterate 0 --nofft input [> output]
$ nwns input [> output]

NW-NS-PartTree-1

3 个参数都设置为最不消耗时间的类型,适合于 ~10,000 到 ~50,000 条序列的比对。

$ mafft --retree 1 --maxiterate 0 --nofft --parttree input [> output]