企业手机网站设计,UE做的比较好的网站,温州专业制作网站,连云港做网站哪里好仅本人练习使用#xff01;#xff01;后续会逐渐修改#xff01;#xff01; mcmctree估算物种分歧时间 - 简书 https://www.cnblogs.com/bio-mary/p/12818888.html 估算系统树分歧时间 —— paml.mcmctree,r8s | 生信技工 http://www.chenlianfu.com/?p2948 4. 使用PAM…仅本人练习使用后续会逐渐修改 mcmctree估算物种分歧时间 - 简书 https://www.cnblogs.com/bio-mary/p/12818888.html 估算系统树分歧时间 —— paml.mcmctree,r8s | 生信技工 http://www.chenlianfu.com/?p2948 4. 使用PAML软件的mcmctree进行物种分歧时间计算
# 新建文件夹并进入
mkdir d.divergence_time cd d.divergence_time
seqfileorthofinder中得到的SpeciesTreeAlignment.fasta文件将物种数和碱基数碱基数将SpeciesTreeAlignment.fasta文件导入TBtools统计每条序列的碱基数即可加在每条序列的前面一行且序列间间隔一行treefileorthofinder中SpeciesTreeAlignment.fasta序列用raxml软件跑出的最优树RAxML_bestTree.jjh_1.tre将后缀改为“trees”并在文件第一行加上物种数目和树的数目先按照上面程序跑一下终端中会出现一个reading master tree然后手动粘贴出来为trees文件在其中加上几个校准点在timetree中查询添加格式参照参考网址。 TimeTree :: The Timescale of Life 另外一定要注明外群的校准点和外群之外的校准点。开头和结尾不用补括号。想看校准点位置加的是否正确可以在Figtree中查看。 mcmctree.ctl路径输入时直接右键复制该文件粘贴在终端里即可
(1) 准备Newick格式的树文件根据RAxML的结果提取系统发育树的拓扑结构
按照paml4.9i/examples中的mcmctree.ctl文件配置自己的数据文件在JJHdata中主要是seq文件和tree文件具体格式参照给出文件。
echo 97 1 input.trees
# 文件内容分两行第一行表述树中有97个物种共计1个树两个数值之间用空分割
# 第二行则是Newick格式树信息其中包含有校准点信息。校准点信息一般指95%HPDHighest Posterior Density对应的置信区间校准点单位是100MYA软件说明文档中使用该单位也推荐使用该单位若使用其它单位后续配置文件中的相关参数也需要对应修改。此外Newick格式的树尾部一定要有分号没有的话程序可能不能正常运行
# echo (((laame,plost),(parub,sccit)0.761.08)2.292.66,(lasul,phgig)); input.trees
# Paxillus rubicundulus Scleroderma citrinum 89MYA 76-108MYA
# Paxillus rubicundulus Laccaria amethystina 246MYA 229-266MYA
# Hymenochaetaceae 125Mya Agaricales-90Mya Agaricomycetes-290Mya Basidiomycota-400Mya
# 把iqtree/raxml-ng等建树软件中得到的树文件中的支持率和枝长信息删除添加化石校准点时间信息格式是时间范围’0.230.26’或者时间点‘0.245’)单位时百万年前100Ma再在首行添加两个数字物种数量和树的数量空格隔开可得到input.tre文件。
手动添加并用Figtree检验
1(JJH) Rickenella Sanghuangporus 227–277Mya
2(JJH) Trametes Heterobasidion 187–236Mya
3(JJH) Pleurotus Schizophyllum 110–201Mya
4(other) Laccaria bicolor Coprinopsis cinerea 120–180Mya
(2) 准备多序列比对Phylip格式输入文件可以综合三个位点核酸序列的多序列比对文件
/media/aa/DATA/JJH/software/standard-RAxML-master/usefulScripts/convertFasta2Phylip.sh ../c.species_tree/allSingleCopyOrthologsAlign.Codon.fasta input.txt
perl -p -i -e s/\s/ / input.txtecho $(cat ../c.species_tree/allSingleCopyOrthologsAlign.Codon.fasta) |sed s/ /\n/g |sed s///g|sed s/ / /g input2.phy
(3) 准备paml mcmctree配置文件
JJH) ndata 1; seqtype 2; burnin 1,000,000; sampfreq 10; nsample 500,000.
# 软件自带测试数据是3组数据而我们的数据只有一组
# 修改mcmctree.ctl
seqfile input.txt
treefile input.trees
mcmcfile mcmc.txt
outfile out.txt
ndata 1 # 输入的多序列比对的数据个数这里密码子3个位置的数据如果有一个则设置为1
seqtype 0 * 0: nucleotides; 1:codons; 2:AAs #数据类型
usedata 3 * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV) # 设置是否利用多序列比对的数据
#0表示不使用多序列比对数据则不会进行likelihood计算虽然能得到mcmc树且计算速度飞快但是其分歧时间结果是有问题的
#1表示使用多序列比对数据进行likelihood计算正常进行MCMC是一般使用的参数;
#2进行正常的approximation likelihood分析此时不需要读取多序列比对数据直接读取当前目录中的in.BV文件。该文件是使用usedata 3参数生成的out.BV文件重命名而来的。
#此外由于程序BUG当设置usedata 2时一定要在改行参数后加 *否则程序报错 Error: file name empty..
#3程序利用多序列比对数据调用baseml/codeml命令对数据进行分析生成out.BV文件。由于mcmctree调用baseml/codeml进行计算的参数设置可能不太好特别时对蛋白序列进行计算时
#推荐自己修该软件自动生成的baseml/codeml配置文件然后再手动运行baseml/codeml命令再整合其结果文件为out.BV文件。
clock 2 * 1: global clock; 2: independent rates; 3: correlated rates
RootAge 8.0 * safe constraint on root age, used if no fossil for root.
model 4
alpha 0.5
ncatG 5
burnin 4000000 #将前x次迭代burnin后再进行取样即打印出该次迭代计算的结果信息各内部节点分歧时间、平均进化速率、sigma2值和各分支进化速率等
sampfreq 10 #每10次迭代则取样一次。
nsample 100000 # 当取样次数达到该次数时则取样结束程序也将运行结束。
(4) 运行PAML软件的mcmctree (若数据量较大则每当程序调用baseml时按ctrl c终止再手动并行化运行baseml)
# 进入conda环境
conda activate FigTree
# 运行
mcmctree mcmctree.ctl
cp out.BV in.BV
(5) 使用mcmctree进行approximate likelihood分析
# 配置文件中参数 usedata 2再次运行同样的mcmctree命令
# 新建文件夹
mkdir run01 run02
cp input.txt input.trees mcmctree.ctl in.BV run01
cp input.txt input.trees mcmctree.ctl in.BV run02
cd run01; mcmctree mcmctree.ctl mcmctree.log
cd run02; mcmctree mcmctree.ctl mcmctree.log (7) 比较两次运行的MCMC树若差异较小则认可其结果
perl -n -e my $out; while (s/(.*?(\d\.\d))//) { my $info $1; my $value $2; my $new $value * 100; $info ~ s/$value/$new/; $out . $info; } $out . $_; print $out run01/FigTree.tre tree01.nex
perl -n -e my $out; while (s/(.*?(\d\.\d))//) { my $info $1; my $value $2; my $new $value * 100; $info ~ s/$value/$new/; $out . $info; } $out . $_; print $out run02/FigTree.tre tree02.nex
perl -e while () { if (s/\s*UTREE.*?\s*//) { s/\s*\[.*?\]//g; print; } } tree01.nex tree01.txt
perl -e while () { if (s/\s*UTREE.*?\s*//) { s/\s*\[.*?\]//g; print; } } tree02.nex tree02.txt
/media/aa/DATA/DATA/SZQ2/command/clf/bin/calculating_branchLength_bias_percentage_of_two_trees.pl --no_normalization_of_total_branch_length tree01.txt tree02.txt bias_of_2runs.txt
# 1.41% 0.77% 0.77%
rm tree*
perl -n -e my $out; while (s/(.*?(\d\.\d))//) { my $info $1; my $value $2; my $new $value * 100; $info ~ s/$value/$new/; $out . $info; } $out . $_; print $out run01/FigTree.tre tree_abbr.mcmctree(8) 将树信息中的简称换成属名和种名全称选做
cp tree_abbr.mcmctree tree_fullName.mcmctree
cut -f 1,2 ../a.preparing_data/source.txt | perl -p -e s/\t/\t\/; s/$/\/; | perl -p -e s#\s#/#; s#^#perl -p -i -e s/#; s#\n$#/ tree_fullName.mcmctree\n#; | perl -p -e s#/\#/\\\#; s#\/#\\\/#; | sh
java -jar /opt/biosoft/FigTree_v1.4.4/lib/figtree.jar tree_fullName.mcmctree
cd ..
最后利用chiplot作图
地质时间标尺网站https://energyeducation.ca/encyclopedia/Geologic_time_scale
