最近碰到了个需求,让我构建某一基因家族的进化树,并根据进化关系进行相应的分类。这让我想起了之前上课的时候,一个做进化的老师给我们讲过,如果不是纯做进化方向的课题,MEGA完全够用了。由于windows的内存等有限,做几十个基因还凑合,要是上百个基因就吃不消了,于是就想到了用Linux下的MEGA来做。
https://www.megasoftware.net/
由于是二进制文件,直接解压缩,添加到环境变量就可以用了。
具体请看我这篇文章。 https://www.jianshu.com/p/5fd60c818651
上一步我得到了该基因家族的所有基因家族的蛋白序列,然后我用windows下的MEGA的 muscle 算法进行了比较,【align-build alignment-上一步的基因家族蛋白序列-muscle比对-data-export-FASTA format】
最终我得到了比对后的 multiproteins.fasta 文件。
首先进行参数的解读,相比于其他软件,我觉得这款软件比较好理解,也容易上手。
这里的 .mao 文件尤为重要,较为简单的方法是拿到windows下去设置,具体请看组学大讲堂的这篇推送。
https://www.omicsclass.com/article/568
版本信息
MEGA version 10.1.8
For 64-bit Linux
Build 10200331
参数解读
EXAMPLES
This example performs a multiple sequence alignment on codons (it assumes that you have created the file "Clustal_Codon_Alignment.mao" using the prototyper (megaproto). A fasta file with coding data is used as input and the resulting alignment is output in the MEGA format:
This example shows how to construct a neighbor-joining phylogeny for each of a list of sequence data files.
The analysis will be performed for each file listed in "listOfDataFiles.txt" and all results will be written to
the ~/Documents/outputDirectory/ directory:
megacc -a ~/Documents/NJ_Tree_Settings.mao -l ~/Documents/listOfDataFiles.txt -o ~/Documents/outputDirectory/
LIST FORMAT
When using the -l option, each file to be analyzed must be on its own line. For example:
~/Documents/myData/seqData1.fas
~/Documents/myData/seqData2.fas
~/Documents/myData/seqData3.fas
If the analyses are to use a user-provided Newick tree file, then the tree files are given on the same line as the data files, following two pipe characters. For example:
~/Documents/myData/seqData1.fas || ~/Documents/myData/treeFile1.nwk
~/Documents/myData/seqData2.fas || ~/Documents/myData/treeFile2.nwk
~/Documents/myData/seqData3.fas || ~/Documents/myData/treeFile3.nwk
我的最终使用:
下一步我打算用 ggtree 来美化,具体学习情况,我再更新。
1. CIPRES网站:
CIPRES Science Gateway V 3.3 : https://www.phylo.org/portal2/home.action#
2. 上传数据,选择RAXML tools
3. 设置参数
首先安装Linux版本RAxML
$ raxmlHPC -f a -x 12345 -p 12345 -# 100 -m PROTGAMMALGX -s example.phy -n ex -T 20
并行化软件支持,能最快速计算。并行化20个任务,每个任务使用8线程,能使用全部160线程计算资源:
$ /opt/sysoft/mpich2-1.5/bin/mpirun -np 20 raxmlHPC -f a -x 12345 -p 12345 -# 100 -m PROTGAMMALGX -s ex.phy -n ex -T 8
参数解读:
RAxML 的参数非常多,设置非常复杂,上述常用例子的参数为:
-f a
此参数用于选择 RAxML 运算的算法。可以设定的值非常之多。 a 表示执行快速 Bootstrap 分析并搜索最佳得分的 ML 树。
-x 12345
指定一个 int 数作为随机种子,以启用快速 Bootstrap 算法。
-p 12345
指定一个随机数作为 parsimony inferences 的种子。
-# 100
指定 bootstrap 的次数。
-m PROTGAMMALGX
指定核苷酸或氨基酸替代模型。PROTGAMMALGX 的解释: "PROT" 表示氨基酸替代模型; GAMMA 表示使用 GAMMA 模型; X 表示使用最大似然法估计碱基频率。
-s ex.phy
指定输入文件。phy 格式的多序列比对结果。软件包中包含一个程序来将 fasta 格式转换为 phy 格式。
-n ex
输出文件的后缀为 .ex 。
-T 20
指定多线程运行的 CPUs 。
RAxML_bootstrap.ex: 完整的全部bootstrapped trees
RAxML_bestTree.ex : 最佳得分 ML 树
RAxML_bipartitions.ex : 有 bootstrap 分值支持的最佳得分树,分值在 node 上。
RAxML_bipartionsBranchLabels.ex : 有 bootstrap 分值支持的最佳得分树, 分值在 branch 上。FigTree不能识别此文件。
说明:
-f a: tell RAxML to conduct a rapid Bootstrap analysis and search for the best-scoring ML tree in one single program run. Example: raxmlHPC -f a -s alg -x 12345 -# 100 -m GTRCAT -n TEST.
f b: when this is specified RAxML will draw the bipartitions using a bunch of topologies (typically bootstrapped trees) specified with -z (see below) onto a single tree topology specified by -t (typically the best-scoring ML tree). Example: raxmlHPC -f b -t ref -z trees -m GTRCAT -s alg -n TEST.
If you use -f a the tree file name RAxML_bipartions.[run_name] will be the best tree with bootstraps. For the f -b approach you first need to estimate a an optimal tree.
Figtree: 使用Figtree查看 RAxML_bipartions.ex 文件,可显示最佳得分树的bootstrap信息。
MEGA: 使用MEGA软件可得到压缩树condensed tree,打开 RAxML_bipartions.ex或者RAxML_bipartionsBranchLabels.ex 显示设定bootstrap阈值的condensed tree。
更多参数可参考:https://www.yuque.com/wusheng/gw7a9p/mcc73y
最近杂事真的非常的满,终于找到时间更新一下。。。。
通过上一篇文章的介绍, 系统发育树的基本概念 大家已经了解清楚,那到底怎么获得一棵可信的进化树呢?
对于群体遗传学分析,一般都会以群体SNPs位点数据构建系统发育树,因此,接下来我主要以SNPs数据为例,介绍系统进化树的构建方法。
序列比对->建树方法选择->计算最佳替代模型->进化树建立->进化树美化
常见的序列比对软件包括:Clustal和Muscle等。
Clustal 除了有自己独立的软件外(多种 *** 作系统都支持),也常被整合到一些常见的软件中,如:Bioedit、MEGA等。
Muscle 同样支持多种 *** 作系统。
两个软件的引用频率都很高,没有绝对的谁好谁坏,哪个顺手就用哪个即可。
1、Distance-based methods 距离法:
基于距离的方法:首先通过各个物种之间的比较,根据一定的假设(进化距离模型)推导得出分类群之间的进化距离,构建一个进化距离矩阵。进化树的构建则是基于这个矩阵中的进化距离关系。
2、Character-based methods 特征法:
基于特征的方法:不计算序列间的距离,而是将序列中有差异的位点作为单独的特征,并根据这些特征来建树。
模型选择的依据如下图:
UPGMA法已经较少使用。一般来讲,如果模型合适,ML的效果较好。对近缘序列,有人喜欢MP,因为用的假设最少。MP一般不用在远缘序列上,这时一般用NJ或ML。对相似度很低的序列,NJ往往会出现Long-branch attraction(LBA,长枝吸引现象),有时严重干扰进化树的构建。贝叶斯方法则太慢。对于各种方法构建分子进化树的准确性,有一篇综述 (Hall BG, 2005) 认为贝叶斯的方法最好,其次是ML,然后是MP。其实如果序列的相似性较高,各种方法都会得到不错的结果,模型间的差别也不大。不过现在文章普遍使用的是NJ是ML模型。
系统发育分析中,最大似然法(ML)和贝叶斯法(BI)是对替代模型非常敏感的两种算法,因此,利用ML法或BI法重建系统发育树前,替代模型的选择是必不可少的过程。
Win *** 作系统下jModeltest的使用方法参考这篇文章: 图解核苷酸替代模型的选择 - jModelTest 篇(By Raindy) 。
ProTest的使用方法可以参考这篇文章: 使用 ProtTest 来选择最优氨基酸替代模型 。
我自己基本都用的是Linux版本的jModelTest,使用及其简单,命令如下:
参数说明:
-d:输入文件。注意!这个软件需要输入的是.phy格式文件,不是.fasta格式。
-f:include models with unequals base frecuencies
-g:include models with rate variation among sites and number of categories
-i: include models with a proportion invariable sites
-s:number of substitution schemes
-v:do model averaging and parameter importances
-a:estimate model-averaged phylogeny for each active criterion
-BIC:calculate the Bayesian Information Criterion
-AIC:calculate the Akaike Information Criterion
结果的最下方,有如图所示的列举,也就是得分最高的模型。
计算完最佳模型,我们就要开始建树了。对于ML树的构建,推荐大家使用新一代RAxML——raxml-ng。
RAxML一直是ML建树的经典工具,其由来自德国海德堡理论科学研究所(Heidelberg Institute for Theoretical Studies)的Alexandros Stamatakis开发。近年来,其江湖地位也受到来自其他软件,尤其是IQ-Tree的挑战。Zhou等人的文章 Evaluating Fast Maximum Likelihood-Based Phylogenetic Programs Using Empirical Phylogenomic Data set 对RAxML,IQ-TREE,FastTree,Phyml四个最大似然法建树软件的实际效果和表现进行了系统比较,其中一个结论是IQTREE在准确性方面要略胜一筹。
近日,RAxML的升级版, raxml-ng 发布!
相较于上一代,raxml-ng有如下优势:
话不多说,直接建树:
参数说明:
--all:Perform an all-in-one analysis (ML tree search + non-parametric bootstrap)
--msa:对其后的序列文件
--model:直接输入上一步产生的最佳模型
--bs-trees:检查树的鲁棒性(robustness)进行自展(bootstrap)检验,进行1000次bootstrapping抽样
--threads:给定线程
运行后结果如下图所示,其中.bestTree就是我们要的树文件,导入树可视化工具即可(我比较常用MEGA和 iTOL ),下次再写一下如何美化进化树吧。
做进化分析的工友们可能都有个感觉,很多分析一等就是好几天,特别是建树(做过的都知道其中的痛苦),有时候忽然加入一个样品又要从头来。因此,一台给力的服务器是必要的工具。比如,上文提到了SNP进化树,我做的还仅仅只是相近物种,而且基因组很小(9M),SNP位点就有4万个,如果要用我MEGA这些软件调用我电脑8核的CPU,1000自展值可能要跑到毕业。
生物学背景出身的我,抄着那一点可怜的计算机常识,在我们课题组购买服务器时,我做了非常多的功课。当然,主要还是听取公司技术人员的建议,通过我非常非常非常长时间的测试,多次使用常见的生物信息分析软件(我主要从事寄生虫基因组、宿主转录组、16S宏基因组等研究),最终,找到了一个性价比超高的服务器配置,具体配置如下:
真心感谢一下烽伟的技术小哥哥们,乐死不疲的回答我各种低级的问题,如果有啥需要可以联系一下他们的技术,感觉蛮靠谱哒,官网: 烽伟科技 。
上一个他们的LOGO,以表感谢。
本文为本人的学习笔记,希望对大家有所帮助。本文大量参考网络文章,文章来源列举于全文末尾。
参考:
一文读懂进化树
使用 ProtTest 来选择最优氨基酸替代模型
RAxML进化树构建的新一代——raxml-ng
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)