步骤一：数据清洗
将ID重复的进行筛除，经过初步计算进化速率为0.9E-3，将同年序列差异<=3的进行随机筛选，保留一个。最后得到147条VP7序列。我为啥做这一步，我在N多文献里就看到个MERS的做了一点这个。其实我想想也可以不做，计算资源够也没关系吧。
步骤二：了解一些基础知识
进化模型：
基于距离有：
• UPGMA（类平均法）
• ME（Minimum Evolution，最小进化法）
• NJ（Neighbor-Joining，邻接法）
基于特征有：
• MP（Maximum parsimony，最大简约法，认为进化历程中发生碱基替代次数最少的系统发育树为最优树）
• ML（Maximum likelihood，最大似然法，基于碱基替代模型，认为似然值最大的系统发育树为最优树）
• Bayesian Inference (BI，贝叶斯方法，在马尔可夫链(Monte Carlo MarkoChain, MCMC)中取样，选取出现频率最高的系统发育树为最优树；BI法通过后验概率(posterior probability, PP)来评估拓扑结构的可靠性，一般PP > 95%的分支可信度较高）
而ML、BI对碱基/氨基酸替代模型敏感，需要对碱基替代模型做评估。通过modelfinder（iqtree内置）、jModeltest2测试，通过似然率（LRT）检验、BIC、AIC、DT指标筛选模型。最终选择GTR+I+G4。
其实各个模型的差异也就那么几个地方：碱基替代频率是否相等、碱基出现的概率是否相等、transition rate (alpha) 、transversion rate (beta)、(PAUP* rate classification）等。有文献表明，模型选择这一步可以不做，直接上GTR+G+I，或者其他模型，最后得出的差异很小的 [Model selection may not be a mandatory step for phylogeny reconstruction]。我可太谢谢这个文章了，因为前期在模型选择上费了好多时间，整理模型之间的差异。所以最后直接选择GTR+I+G和HKY+I+G。因为看到很多文献也蛮喜欢用这2个的，还有用JC的。但没有太多文献具体说明为啥选择这个模型，最多就是说经过jmodeltest筛出BIC最小的。
步骤三：配置文件
BEAST是用的比较多的一个软件，先要配置好xml文件，其实就是模型标记文件（我记得python的pickle包存模型也是这个格式，于是我测试下是否pickle也识别，哈哈哈，当然不能啦！笨笨笨！）。BEAST是一个软件集，先用BEAUTi进行配置。将比对好的fa文件拖入：
Tips选项，可以把时间自动识别。因为我的ID最末是时间，根据自己的格式处理：
BEATUi tips
Sites选项：这个其实就是位点的意思，每一个碱基/氨基酸：替代模型我这选择的GTR，base freq选择est、G4、codon那选择的(1+2),3 。因为密码子的第三个碱基相对前2个变化要大。可以设置为off、1,2,3看看xml文件里是怎么描述模型的：
codon1 2 3的变异速率不一致
BEAUTi Site
off
（1+2），3
1，2，3[比前者多一个3到3]
Trees选项：种群增长模型，我选择的贝叶斯天际线BS
clocks 选项：看你想用什么分子钟，我选择的严格分子钟。如何判断自己的数据适合宽松的分子钟还是严格的分子钟？（此处复制From高芳銮教授QQ空间日志教程）
方法一：可以先用uncorrelated relaxed lognormal clock配合默认参数配置个xml文件，MCMC链长不需要很大，运行结束后，将xx.log文件通过Tracer导入查看coefficientOfVariation（或ucld.stdev）值，如果该值接近0，说明数据适合严格分子钟，反之说明数据适合宽松分子钟。
方法二：通过贝叶斯因子法（Bayes Factor, BF）来比较分子钟模型和树先验（tree prior）模型之间的组合的边际似然值(Marginal Likelihood)，通过BF法分析后，哪个模型对应的LnL值越大，相对更优,此方法既可以确定分子钟模型，也可以用于确定哪种树先验模型更加适合自己的数据。具体也可以通过以下操作来确定：在Tracer 中同时导入两个模型的xx.log 文件，选中所有要比较的模型后，依次点击菜单Analysis--->Model Comparation，Analysistype改为Harmonic mean，Reps=100，最后根据获得的BF值（BF=Pr(D|M1)/Pr(D|M2)至少大于3，则说明M1>M2）范围判断模型间的优劣程度。
方法三（更新）: 同样是采用比较分子钟模型和树先验（tree prior）模型之间的组合的边际似然值(Marginal Likelihood)，但是不使用BF法，而是采用path sampling（路径采样）分析，相比BF法，path sampling更加精准。（参考文献：Improving the accuracy of demographic and molecular clock model comparison while accommodating phylogenetic uncertainty)。
priors 选项：我理解的就是对模型初始的分布函数进行设置，比如函数形状、mean、sd等。主要是对你前面进行设置的做模型先验函数预估。
MCMC选项：length of chain ，看资料说是不小于1000*序列数**2。我是147条，我直接设置的200000000，跑了我一天，ESS全部>200，开森！因为前面有设置小一点，有group的ESS<200。所以遇到有的参数不收敛的时候，可以尝试把chain长度设置大一点！另外的参数自己看字面也能理解，不多说。
MCMC
r然后生成xml文件。再丢到服务器上，跑beast。
one thousand years later......
得到了以下4个文件：待tree是树文件啦，自己找个软件/包去美化下，工具很多，在线的也有：
beast结果
把log文件拖到tracer里看看：可以看到进化速率mean值是1.0021E-3/site/year。跟我看的文献里的差不多。
当然，除了BEAST还有iqtree、MEGA都可以做。文献里基本是这几个。iqtree采用极大似然法构建进化树，很快，但我为了准确，-b 1000，然后从昨天跑到我现在写这个文章还没跑完。有manual，可以看看，反正我的命令是：/project/baowenjuan/APP/iqtree/iqtree-1.6.12-Linux/bin/iqtree -s merged_seq_y.mafft.fa -pre iqtree_K3PuFG4_RV -b 1000 -nt AUTO -m K3Pu+F+G4 （这里的模型是iqtree前期评估的最优模型）
jModelTest对模型评估以及构建进化树，自行看manual，都不难。我的命令是：java -jar /project/baowenjuan/APP/jmodeltest-2.1.10/jModelTest.jar -d merged_seq_ymafft.nex -f -g 4 -AIC -AICc -BIC -DT -i -s 11 -a -w -p >jModelTest.o
结果里有它自己根据BIC值选择的模型构建的进化树。
一个礼拜，从0走到现在，一套流程走完啦。具体发现了什么，就需要后续的分析啦！
PS. 今年流年不利，几个项目被砍。肿瘤玩厌了，又搞早筛，现在又“被迫”琢磨了病毒进化。所以接下来我要玩啥呢？逼着我变全能吗？我不想啊，我想专一一点，我想往大数据、建模方向发展呢！加油！
网友评论